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Preface 


This  study  began  as  an  attempt  to  compare  a  new  method 
for  predicting  stability  with  the  coefficient  method.  The  new 
method,  using  Pade ’  approximants  and  systems  of  linear 
equations,  had  predicted  the  oscillatory  behavior  discovered 
by  a  previous  thesis  (8;  14:246).  To  do  a  comparative  study, 
a  new  approach  using  the  space  and  time  step  domain  was 
developed.  Also,  for  a  comparative  study  it  was  felt  that  a 
survey  of  many  finite  difference  methods  for  both  the  one  and 
two  dimensional  cases  was  required,  so  this  thesis  looks  at 
five  popular  finite  difference  methods. 

I  would  like  to  thank  Dr.  Bernard  Kaplan  of  the  Air  Force 
Institute  of  Technology  for  letting  me  pursue  this  independent 
study  as  just  that,  an  independent  study,  while  always 
demanding  quality  results.  His  guidance  kept  me  from  getting 
too  far  off  track  at  times.  Also,  I  would  like  to  express  my 
thanks  to  Dr.  N.  Pagano  of  AFWAL/MLBC  for  sponsoring  this 
study.  Finally,  I  would  like  to  express  my  thanks  to  Major 
Lupo,  whose  door  was  always  open,  especially  when  it  came  to 
questions  about  mainframe  computers,  and  without  whose  support 
the  graphics  in  this  thesis  would  still  just  be  numbers 
floating  around  somewhere  in  the  VAX  computer. 


i  For 


—  Thomas  M.  Dipp 
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DF  Maximum  Numeric  Time  Oscillations,  IM  Start  . 
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Figure 

21.  CN  Maximum  Numeric  Time  Oscillations 
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Notation  and  Abbreviations 


=  Ay  =  h  =  spatial  increment  between  nodal  points. 

s  Ny  =  N  =  total  spatial  node  points  in  x  or  y  direction. 

L  =  Nh  =  total  length  of  rod  or  square  plate  in  x  or  y 
direction. 

At  *  k  *  time  increment  (seconds) 

t  =  time  (seconds) 

M  =  total  time  nodes. 

T  =  Mk  =  total  time  to  solution 

a  =  thermal  diffusivity,  here  always  =  1.0 
r  =  a£t/Ax2  =  aAt/Ay2  =  ak/h2 

U(i,j,n)  =  temperature  at  the  ith,  jth  node  at  time  step  n, 
=  U(iAx,  jAy.nAt)  =  U(x,y,t) 

[B]  =  notation  for  matrix  B 

<b>  =  notation  for  column  matrix  b 

EX  =  Fully  Explicit  Method 

DF  =  DuFort-Franke 1  Method 

CN  =  Crank-Nicol son  Method 

ADI  =  Peaceman-Rachf ord  Alternating  Direction  Implicit 
Method 

IM  =  Fully  Implicit  Method 

1- D  =  one  dimensional 

2- D  =  two  dimensional 

FDM  =  Finite  Difference  Method 

[A]  *  Pade ’  exponential  matrix  A 

[G]  =  amplification  or  growth  matrix  G 


c 


s=l(l)n  =  s  varies  from  1  to  n  by  increments  of  (1),  i.e. 
s  ~  if  2 1  3  f  • • • f  n  if  n 

Zs .  Z« , n  =  1-D  and  2-D  matrix  A  eigenvalues 

Gs >  Ga , n  =  1-D  and  2-D  amplification  matrix  eigenvalues 


=  3.141592654 . . . 


Abstract 


This  thesis  compared  the  coefficient  method  with  the  full 
matrix  method  for  predicting  stability  and  oscillatory 
behavior  of  Finite  Difference  Methods,  FDMs ,  used  in  solving 
the  one  and  two  dimensional  transient  heat  (diffusion) 
equation  with  Dirichlet  boundary  conditions.  Five  FDMs  were 
used:  the  fully  implicit,  fully  explicit,  DuFort-Franke 1 , 

Crank-N i co 1  son ,  and  Peaceman-Rachf ord  ADI. 

Analytically,  the  Pade ’  method  was  found  to  be  equivalent 
to  the  matrix  method  in  predicting  stability  and  oscillations. 
The  matrix  method  was  shown  to  be  more  severe  than  the 
coefficient  method  in  predicting  both  stable  and 
non-osci 1 1 atory  step  size  constraints.  Also,  the  matrix 
method,  while  mathematically  more  rigorous,  proved  to  be  more 
difficult  to  derive  and  analyze,  possibly  limiting  its 
usefulness.  Ultimately,  it  was  found  that  the  coefficient 
method’s  predictions  could  be  derived  from  the  matrix  method 
for  two-level  FDMs. 

Numerically,  all  methods  were  solved  repeatedly  and  the 
results  investigated  for  oscillatory  behavior  and  maximum 


I  errors  in  the  h-k,  or  space  and  time  step,  domain.  These 

t 

c 

‘  results  allowed  functional  comparison  with  stability  and 

oscillatory  predictions  made  by  the  coefficient  and  full 

|  matrix  methods.  Neither  the  coefficient  nor  matrix  method’s 

P 

L*  stability  and  oscillatory  predictions  were  significantly 

t" 
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violated  for  two-level  FDMs .  Finally,  the  matrix  method  could 
generally  predict  when  non-d i sastrous  oscillations  would 
occur,  based  on  moderate  step  size  constraints. 


PREDICTING  OSCILLATORY  FINITE  DIFFERENCE  SOLUTIONS 
TO  THE  HEAT  EQUATION:  A  COMPARATIVE  STUDY 
OF  THE  COEFFICIENT  AND  MATRIX  METHODS 

I .  Introduction 


Background 

Finite  difference  methods  (FDMs)  have  become  increasingly 
important  to  the  advancement  of  science  and  technology  as  ever 
faster  and  more  powerful  digital  computers  are  developed  to 
handle  the  complex  and  often  analytically  intractable  problems 
arising  in  science  and  industry.  One  such  problem,  the 
transient  heat  (diffusion)  partial  differential  equation,  is 
often  solved  using  finite  difference  methods.  However,  as 
with  all  numerical  methods,  the  behavior  of  a  finite 
difference  method  used  to  numerically  solve  a  problem  must  be 
understood  before  the  method  can  be  applied  with  the 
reasonable  hope  of  getting  anything  more  than  useless 
significant  digits  from  a  computer.  In  the  relatively  recent 
past,  methods  proven  to  be  unconditionally  stable  have  been 
embraced  only  to  find  that  unacceptable  errors  and 
oscillations  were  occurring,  such  as  in  the  Crank— N i co 1  son , 
DuFort-Franke 1 ,  and  Peaceman-Rachf ord  alternating  direction 
implicit  FDMs. 


Stability  is  a  condition  where  errors  introduced  into  the 
FDM  are  bounded  as  time  progresses.  Von  Neuman  and  Goldstein 
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note  four  major  sources  of  errors  in  describing  physical 
systems:  1)  Mathematical  models  of  nature  are  approximations, 

2)  Most  descriptions  require  measured  data,  3)  truncation 
error  due  to  most  descriptions  involving  an  infinite 
continuum,  and  limiting  processes  cannot  be  completed,  and  4) 
roundoff  error  caused  by  the  fact  that  calculations  cannot  be 
carried  out  with  an  infinite  number  of  digits  (15:6).  The 
last  source  of  error  listed  above  is  the  most  important  source 
in  the  study  of  FDM  stability  and  results  from  the  fact  that 
finite  difference  methods  involve  limited  precision  in 
computer  arithmetic.  Consequently,  errors  are  introduced  at 
each  time  step  of  the  calculations,  and  would  grow  without 
bound  if  the  FDM  were  not  stable.  Even  if  a  method  is  stable, 
errors  and  oscillations  that  would  eventually  damp  out  with 
time,  may  produce  unacceptable  inaccuracies  in  a  given 
solution.  As  a  result,  a  knowledge  of  techniques  for 
predicting  stability  and  oscillations,  and  how  well  these 
techniques  work,  is  needed  before  solving  a  problem  on  the 
computer . 

The  Two  Dimensional  Transient  Heat  (Diffusion)  Equation 

The  heat  equation  is  a  parabolic  partial  differential 
equation,  and  in  two  dimensions  has  the  form 

b ( x , y , t ) t  =  a  [  U(x,y,t)xx  +  U(x,y,t)yy  ] 


o 


where  U(x,y,t)  is  the  temperature,  ’a’  is  the  thermal 
diffusivity,  t  is  the  time,  and  x  and  y  are  the  spatial 
dimensions  of  length.  Partial  differentiation  of  U(x,y,t) 
with  respect  to  t,  x,  and  y  is  denoted,  in  standard  notation, 
by  subscripted  variables. 


Problem  Statement 

The  goal  of  this  thesis  is  to  compare  the  coefficient  and 
full  matrix  (Pade’  approximant)  methods  of  predicting 
stability  and  oscillatory  behavior  of  FDMs  used  in  solving  the 
transient  heat  (diffusion)  equation. 


Scope 

Five  finite  difference  approximations  to  the  transient 
heat  (diffusion)  equation  are  used  in  this  study: 


Fully  Implicit  (IM) 

Fully  Explicit  (EX) 

DuFort-Franke 1  (DF) 

Crank-Nico 1  son  (CN) 

Peaceman-Rachf ord  Alternating  Direction 
Implicit  (ADI) 


Both  the  one  (i-D)  and  two  (2-D)  dimensional  FDMs  are 
analyzed,  except  in  the  case  of  the  ADI  method  which  has  no 
1-D  representation.  Also,  the  one  and  two  dimensional  FDMs 
are  each  used  to  solve  a  different  test  problem.  Both 


problems  are  boundary  value  problems  using  Dirichlet  boundary 
conditions,  both  assume  a  thermal  diffusivity  of  one,  and  both 
the  rod  and  the  square  have  unit  lengths.  For  the  two 
dimensional  case,  only  the  standard  Ax=Ay  plane  is  used.  The 
number  of  spatial  nodes  varies  from  three  to  25  and  the  number 
of  time  steps  varies  from  four  to  100.  All  nodes  and  all  time 
steps  are  included  in  the  oscillation  and  error  analysis. 

Approach 

First,  two  previous  theses,  where  oscillations  were  found 
to  exist  in  the  ADI  (8:56)  and  DF  (4:62)  FDMs ,  were  reviewed. 
Next,  published  oscillatory  solutions  were  reproduced  for  1-D 
cases  in  Smith  (13:124),  Twizell  (14:206),  and  Lawson  and 
Morris  (10:1217).  From  this  experience,  it  was  realized  that 
no  h-k  domain  studies  comparing  stability  and  oscillation 
predictions  had  been  done.  A  literary  search  was  conducted 
and  no  comparative  domain  studies  were  found.  Two  problems 
with  known  analytical  solutions,  which  were  previously  used  in 
showing  oscillatory  behavior,  were  selected  for  indepth 
stability  analysis  in  the  h-k  domain  ( 1 3 : 1 24 ; 4 : 1 9 ; 8 : 19 ) .  For 
such  a  comparison  study,  a  survey  of  many  FDMs  in  the  h-k 
domain  was  indicated.  Computer  modules  to  determine  stability 
components  were  developed  and  computer  programs  for  the  one 
and  two  dimensional  FDMs  were  written.  Numerical  solutions  to 


the  problems  using  the  FDMs  were  graphed  in  the  h-k  domain  and 
compared  with  predicted  functional  restrictions  on  h  and  k 
obtained  from  the  coefficient  method  and  Pade ’  approximant 


(full  matrix)  method.  Finally,  the  analytic  development  of 
the  coefficient,  Pade '  approximant,  and  full  matrix  methods 
were  explored  to  determine  the  relationships  among  the 
methods . 

Sequence  of  Presentation 

Chapter  2  introduces  general  finite  difference 
approximations  and  the  finite  difference  equations  for  the 
five  FDMs  used  in  this  thesis.  Stability,  oscillations,  and 
accuracy  are  discussed  in  Chapter  3,  along  with  the 
coefficient  method  and  full  matrix  (Pade’  approximant)  method 
for  predicting  stability  and  oscillations.  Because  of  the  new 
approach  to  graphing  FDM  stability  components  and  errors  in 
the  h-k  domain,  which  uses  2  dimensional  contour  plots  of 
three  dimensions,  Chapter  4  covers  the  rationale  behind  the 
domain  approach,  what  it  includes,  how  to  read  the  domain 
graphs,  and  comparisons  with  conventional  stability  graphs. 
Chapter  5  presents  the  numerical  results  for  maximum  time  and 
spatial  oscillations,  as  well  as  for  maximum  absolute  errors, 
as  well  as  comparisons  with  stability  predictions.  Finally, 
conclusions  and  observations  are  given  in  Chapter  6. 

Computer 

A  VAX  11/780  32  bit  computer  was  used.  Programs  were 


written  in  FORTRAN  77.  All  calculations  were  done  using 
double  precision  to  reduce  the  effects  of  round  off  error. 


Extensive  use  of  structured  programming  technique  was  applied 
to  insure  that  code  was  readable  and  modular,  allowing  quick 
program  design  with  reliability.  Graphs  were  drawn  using 
standard  library  routines  provided  by  Major  Lupo  (11),  and 
hardcopy  output  was  obtained  with  an  Imagen  Laser  Printer. 


Finite  Dif ference  Approximations 


All  finite  difference  methods  involve  approximations  to 
differential  equations.  Since  the  analytic  solution  for  the 
heat  equation  for  all  space  and  time  may  be  impossible  or  not 
practical  to  derive,  numerical  solutions  based  on  a  knowledge 
of  the  present  state  of  a  system  and  the  rates  of  change  to 
that  system  must  be  developed  instead.  These  numerical 
solutions  are  accomplished  by  using  approximations  to  the 
partial  derivatives  in  the  governing  partial  differential 
equation,  and  involve  replacing  the  actual  derivatives  by 
approximate  derivatives  built  from  truncated  Taylor  series 
expansions.  The  details  of  this  process  can  be  found  in 
numerous  texts  (1:532-538,578-597;  13:6-24;  14:31-37). 

A  fundamental  result  is  that  space  and  time  are  no  longer 
an  infinite  continuum,  but  become  discretized  by  the 
truncation  of  the  derivative  expressions.  In  effect,  the 
result  is  a  finite  grid  of  difference  approximations  based  on 
adjacent  time  and  space  point  values,  or  nodes.  Inevitably, 
it  is  this  discretization  that  introduces  truncation  error 
into  the  numerical  solution. 

Five  finite  difference  approximations  commonly  used  to 
solve  the  heat  equation  are  the  fully  explicit,  fully 
implicit,  Crank-Nicolson,  DuFort-Franke 1 ,  and 

Peaceman-Rachf ord  alternating  direction  implicit  methods.  In 


c 


listing  the  finite  difference  approximations,  the  following 
notation  is  implemented: 


C 


C 


AX  «  Ay  =  h ,  At  =  k ,  x  =  ih,  y=jh,  t  *nk  where 
i ,  j ,  and  n  are  integers  and 
U (x ,y , t )  =  U ( ih , jh , nk)  =  U(i,j,n)  with 
a  =  1 ,  r  =  ak/h2 . 

Fully  Explicit  (EX)  Approximation 

One  of  the  easiest  methods  to  implement,  the  EX 
calculates  the  temperatures  at  each  new  time  step  only  in 
terms  of  the  old  temperatures  at  the  previous  time  step.  The 
formulas  for  one  and  two  dimensions  are  (9:155,235): 

U ( i , n+ 1 )  *  (l-2r)U(i,n)  +  r[U(i+l,n)  +  U(i-l,n)]  (1) 

U ( i , j , n+1 )  =  ( l-4r)U( i , j ,n)  +  r[U(i+l,j,n)  (2) 

+  U(i-l,j,n)  +  U ( i , j  + 1 . n )  +  U ( i , j - 1 , n ) ] 


DuFort-Franke 1  (DF)  Approximation 

The  DuFort-Franke 1  approximation  is  a  three  level 
explicit  method  that  overcomes  the  instability  of  the  fully 
explicit  method  (9:238).  It  was  the  goal  of  the  DF  method  to 
have  the  ca 1 cu 1  at i ona 1  ease  of  an  explicit  method  without  the 
instability  associated  with  the  fully  explicit  method.  Forms 
for  the  one  and  two  dimensional  approximations  are 
(9: 157,239) : 


N^.VfV.'W. 
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( l+2r)U( i ,n+l )  »  2r[U(i+l,n)  +  U(i-l,n)]  (3) 

+  (l-2r)U(i,n-l) 

( l+4r)U( i , j ,n+l )  *  2r[U(i+l , j ,n)  +  U(i-l,j,n)  (4) 

+  U(i,j+l,n)  +  U(i,j-l,n)J  +  ( l-4r )U( i , j ,n-l 

Crank-Nicol son  (CN)  Implicit  Approximation 

The  Crank-Nicolson  method  equally  weights  and  combines 
the  fully  explicit  and  fully  implicit  methods.  The  goal  of 
the  CN  method  was  to  reduce  the  truncation  error  of  the  IM 
method  and  hopefully  improve  accuracy.  The  one  and  two 
dimensional  forms  are  (9:160,242): 


( 2+2r )U( i , n+1 )  -  r[U(i+l,n+l)  +  U(i-l,n+l)]  (5) 

*  ( 2-2r )U( i , n)  +  r[U(i+l,n)  +  U(i-l,n)] 

(  2+4r  )li (  i  ,  j  , n+1 )  -  r  [U (  i  +  1 ,  j  , n+1 )  +  U(  i-1 ,  j  , n+1 )  (6) 

+  U(i,j+l,n+l)  +  U ( i , j-1 , n+1 ) ] 

=  (2-4r)U(i , j ,n)  +  r[U(i+l,j,n)  +  U(i-l,j,n) 

+  U(i,j+l,n)  +  U(i,j-l,n) 

Peaceman-Rachf ord  (ADI)  Approximation 

The  ADI  method  also  combines  the  EX  and  IM  methods, 
similar  to  CN  method.  However,  the  ADI  method  alternates  in 
its  use  of  the  EX  and  IM  methods.  ADI  solves  a  given  time 
step  explicitly  in  one  direction  and  implicitly  in  the  other, 


9 


and  then  reverses  the  EX  and  1M  directions  in  solving  for  the 


next  time  step.  Both  time  steps  must  have  the  same  duration 


and  both  should  be  completed  for  a  stable  solution.  The  goal 


of  the  ADI  method  was  to  have  the  speed  of  tridiagonal 


matrices  instead  of  the  time  problems  of  the  IM  method.  It 


has  the  following  form  in  two  dimensions  (9:246-247): 


( I  +  2r )U ( i , j , n+1 )  -  r [U( i  +  1 , j , n+1 )  +  U( i-1 , j , n+1 ) ]  (7a) 


=  ( l-2r )  U(i,j,n)  +  r[U(i,j+l,n)  +  U(i,j-l,n)] 


to  advance  from  time  step  n  to  n+1,  followed  by 


( l  +  2r  )li  (  i  ,  j  ,  n+2  )  -  r  [U  (  i  ,  j  +  1  ,  n+2  )  +  U (  i  ,  j-1  ,  n+2  ) )  (7b) 


=  ( l-2r)U( i , j ,n+l )  +  r  [U ( i  +  1 , j , n+1 )  +  U ( i-1 , j  ,  n+1 ) ] 


to  advance  from  time  step  n+1  to  n+2. 


Fully  Implicit  (IM)  Approximation 


By  replacing  the  second  order  spatial  derivatives  with 


the  second  order  backward-difference  approximation,  this 


method  computes  the  temperatures  at  the  new  time  step  based  on 


temperatures  at  neighboring  node  points  which  are  also  at  the 


new  time  step.  The  one  and  two  dimensional  form  are 


(9: 159,242) : 


( l+2r )U( i , n+1 )  -  r[u(i+l,n+l)  +  U(i-l,n+l)]  =  U(i,n)  (8) 


( l+4r )U( i , j , n+1 )  -  r [U( i+1 , j , n+1 )  +  U( i-1 , j , n+1 ) 
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+  U(  i  ,  j  +  1 , n+1 )  +  U ( i , j -1 , n+1 ) ]  =  U(i,j,n) 

Methods  of  Solution 

The  EX  and  DF  methods  are  explicit  and  are  readily  solved 
with  one  equation  for  each  unknown.  However,  the  DF  method 
requires  two  initial  time  steps  be  known  before  it  becomes 
explicit.  Therefore,  the  DF  method  is  often  started  by  using 
another  method  to  calculate  the  first  time  step  from  the 
initial  and  boundary  conditions.  Initial  time  steps  were 
calculated  using  the  analytic  solution  and  the  IM  method, 
because  they  do  not  introduce  oscillations  into  the  first  time 
step.  However,  the  one  dimensional  case  was  also  started 
using  the  CN  method  for  comparison  with  an  FDM  that  could 
osci 1 1  ate . 

The  C.\ ,  ADI,  and  IM  are  all  implicit  formulations,  and 
require  that  a  system  of  simultaneous  equations  be  solved  in 
order  to  calculate  the  unknowns.  The  systems  of  equations 
have  the  general  form: 

[A] < x >  *  <b> 

where  [A]  is  an  NxN  matrix  for  the  one  dimensional  case  and  an 
N2xN2  matrix  for  the  two  dimensional  case,  <x>  is  a  column 
vector  composed  of  the  new  temperatures  at  t+k,  <b>  is  a 
column  vector  containing  the  old  temperatures  and  boundary 
conditions,  and  N  =  L/h  ,  where  L  =  X2-X1.  xi  and  X2 


1 1 


represent  the  length  of  a  1-D  rod  or  the  side  of  a  2-D  square 
plate. 

One  dimensional  solutions  of  the  column  vector  <x>  for 
the  CN  and  IM  methods  were  obtained  efficiently  and  accurately 
by  using  a  tridiagonal  solver.  Since  the  ADI  method 
alternately  solves  using  the  EX  and  IM  methods,  the  ADI  method 
also  obtained  solutions  for  <x>  by  using  a  tridiagonal  solver. 
Finally,  the  two  dimensional  CN  and  IM  methods  were  solved  to 
at  least  13  significant  digits  of  accuracy  by  using 
Gauss-Seidel  Iteration  combined  with  Successive 
Over-Relaxation  (1:428,434). 

At  least  10  digits  of  accuracy  and  precision  was 
constantly  maintained  in  the  numerical  solutions  of  the  FDMs 
since  1)  double  precision  was  used  exclusively  in  all  the 
finite  difference  programs,  2)  matrix  sizes  for  direct  matrix 
methods  did  not  exceed  25x25  with  no  more  than  100  time  steps, 
and  3)  two  dimensional  implicit  matrices  were  solved  using 
Gauss-Seidel  iteration  until  13  digits  of  numerical 
convergence  were  achieved. 
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III.  Stability,  Osci I lations ,  and  Accuracy 


Stability  is  a  property  of  FDMs  which  guarantees  that  the 
numerical  solution  to  the  finite  difference  equation  will  not 
grow  without  bound  as  time  progresses.  As  such,  stability  is 
a  property  that  is  independent  of  the  actual  partial 
differential  equation.  Since  numerical  error  introduced 
during  the  calculations,  mainly  from  roundoff  error,  can  be 
treated  as  perturbations  to  the  vector  of  initial  values, 
numerical  error  is  subject  to  the  same  arithmetic  operations 
as  the  vector  of  initial  values  (13:48,53).  Therefore, 
stability  also  guarantees  that  numerical  error  is  bounded 
(3:98) . 

Stability  has  also  been  given  more  specific  definitions, 
such  as  conditional  or  unconditional  stability,  and 
Ao-stability  or  Lo-stabi  1  i ty .  Conditional  stability  refers  «.o 
a  finite  difference  equation  that  is  only  stable  for  certain 
values  of  step  sizes  for  time  and  space,  and  unstable  for  all 
other  step  sizes  (9:167).  Unconditionally  stable  FDMs  have  no 
step  size  restrictions,  and  are  stable  over  the  whole  h-k 
domain.  Ao-stability  guarantees  unconditional  stability  while 
Lo-stability  not  only  guarantees  unconditional  stability,  but 
also  guarantees  no  oscillations  as  well  (14:213-215).  These 
definitions,  when  applied  to  the  five  types  of  FDMs  use  in 
this  thesis,  are  as  follows: 


EX 

Conditional  ( 1 

:581),  Neither  Ao  nor  Lo 

(13:122) 

DF 

Unconditional 

(13:153),  Ao 

(definition, 

13:120) 

CN 

Unconditiona 1 

(13:56),  Ao 

(13:121) 

ADI 

Unconditional 

(definition, 

13:120),  Ao 

(14:240) 

IM 

Unconditional 

(1:582),  Lo 

(13:121) 

However,  stability  does  not  guarantee  accuracy.  This  is 
because  stability  is  a  necessary,  but  not  sufficient, 
condition  for  accuracy  (9:167).  Also,  a  FDM  may  be  stable, 
yet  still  experience  oscillations  which  can  have  disastrous 
effects  on  the  accuracy  of  a  solution  (12:281-84).  Further, 
the  stability  of  an  FDM  can  depend  on  other  factors,  such  as 
the  type  of  boundary  conditions  employed  (9:167).  Still, 
because  stability  is  necessary  to  accuracy,  it  is  an  important 
first  step  in  analyzing  any  FDM.  Accuracy  and  computational 
efficiency  must  become  secondary  considerations  (9:167). 

Two  methods  of  analyzing  stability  and  oscillations  were 
investigated  and  used:  the  coefficient  method  and  the  Pade ’ 
(full  matrix)  method.  Their  description  follows,  beginning 
with  the  coefficient  method. 

Coefficient  Method 

The  coefficient  method  is  an  easy  to  use  and  understand 
method  of  predicting  stability  and  oscillations  in  FDMs .  In 
this  method,  the  FDM  behavior  for  all  of  the  nodal  points  is 
approximated  by  applying  the  finite  difference  equations  to  a 
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one  nodal  point  grid  with  homogeneous  boundary  conditions. 
Consequently,  since  only  one  nodal  point  is  examined,  the 
coefficient  method  is  often  considered  a  crude  one  nodal  point 
approximation  when  applied  to  predicting  the  stability  of  all 
the  nodal  points. 

Myers,  in  his  book,  Analytical  Methods  in  Conduction  Heat 
Transfer ,  provides  a  clear  and  concise  look  at  the  one 
dimensional  case  of  the  coefficient  method  (12:282).  This 
method  has  been  directly  extended  to  the  two  dimension  heat 
equations  in  previous  papers  on  the  subject  of  stability  and 
oscillations  (8:10-14,4:10-12).  Resulting  ratios  for  the  one 
and  two  dimensional  cases,  Equations  (1)  through  (9),  are 
listed  below: 


Ful ly  Expl icit 

1-D : 

P 

= 

l-2r 

Ful ly  Expl icit 

2-D: 

P 

= 

l-4r 

DuFort-Franke 1 

1-D: 

P 

= 

(l-2r)/(l+2r) 

DuFort-Franke 1 

2-D: 

P 

= 

( l-4r ) / ( l+4r ) 

Crank-N i co 1  son 

1-D: 

P 

= 

(l-r)/(l+r) 

Crank-Nicolson 

2-D: 

P 

= 

( 1 -2r ) / ( 1 +2r ) 

Peaceman-Rachf ord  2-D: 

P 

= 

( l-2r ) / ( l+2r) 

Fully  Impl icit 

1-D: 

P 

= 

l/( l+2r) 

Fully  Implicit 

2-D: 

P 

= 

l/( l+4r) 

where  p  is  the  ratio  of  the  new  to  old  temperatures  at  the 
single  node  point,  and  is  defined  as: 
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l-D :  p  =  U ( i , n+1 ) /U ( i , n ) 


2-D:  p  =  U( i , j , n+1 )/U( i , j , n) 


Stability  and  oscillatory  behavior  of  a  finite  difference 


method  can  then  be  inferred  from  the  various  values  the 


p-ratio  takes  on  as  the  value  of  r  changes,  where  r=k/h2  is  a 


positive  number  ranging  from  zero  to  infinity.  This  behavior 


is  described  bj-  Myers  and  in  summary  is  as  follows  (12:282): 


p  >  1  Unstable,  No  oscillations  =  unstable  growth 


0  <  p  <  1  Stable,  No  oscillations  =  stable  decay 


-1  <  p  <  0 


Stable,  Oscillations  =  stable  oscillations 


p  <-l 


Unstable,  Oscillations  *  unstable 
osci 1 1  at  ions 


Values  of  the  p-ratio,  and  hence  the  predicted  behavior 


of  Equations  (1)  through  (9),  can  be  calculated  algebraically 


or  graphically  for  various  values  of  r.  Figure  3  (top)  is  a 


graph  that  has  been  modified  to  show  the  behavior  of  all  the 


previously  listed  p-ratios  as  a  function  of  r  (8:14).  From 


the  graph,  it  is  easily  seen  that  the  EX  method  is  the  only 


predicted  conditionally  stable  method,  that  the  IM  method  is 


the  only  predicted  Lo-stable  method,  and  that  the  rest  are 


unconditionally  stable  but  can  display  conditional  oscillatory 


behavior  depending  on  the  value  of  r.  Figure  1  (top) 


summarizes  the  predicted  behavior  by  using  a  bar  chart  to 
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Figure  1.  Coefficient  and  Matrix  Method 
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display  all  the  various  behaviors  versus  common  values  of  r 
where  the  behavior  of  the  FDMs  change.  The  values  of  r  where 
changes  are  predicted  for  various  of  the  FDMs  are  r=l/4, 
r=l/2,  and  r=l.  Finally,  notice  that  the  1-D  DF ,  2-D  ADI,  and 
2-D  CN  all  share  the  exact  same  p-ratios,  and  hence,  same 
predicted  behaviors. 
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Full  Matrix  Method 

The  full  matrix  method  used  to  predict  stability  and 
oscillatory  behavior  is  similar  to  the  coefficient  method,  but 
uses  the  entire  matrix  to  determine  the  behavioral 
characteristics  of  the  finite  difference  equation.  As  such, 
it  is  much  more  complicated  and  requires  a  knowledge  of 
matrices,  matrix  eigenvalues  and  eigenvectors,  and  matrix 
stability.  There  are  two  basic  ways  to  approach  matrix 
stability  which  are  equivalent  but  involve  generating  the 
amplification  matrix  differently.  The  most  straightforward 
approach  is  to  generate,  in  matrix  form,  the  simultaneous 
system  of  equations  necessary  to  describe  the  FDM.  For  a 
simple  two-level  FDM,  where  two  time  levels  are  involved  in 
the  numerical  calculations,  the  resulting  matrix  that  is 
generated  is  the  amplification  matrix.  The  amplification 
matrix,  [G] ,  must  be  in  the  form  where 


<U(i,j,n+l)>  =  [G] <U(i , j ,n)> 

which  can  be  easily  and  directly  achieved  for  explicit 
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methods,  but  requires  for  implicit  forms,  such  as  CN,  the 
following  set  of  manipulations: 


[B]<Un+i>  =  [C]<l!n>  be  changed  to 
<U(i,j,n+l)>  =  [B] - » [C] <U( i , j , n) >  or 
<U ( i , j , n+1 )  >  =  [G)<U(i,j,n)> 

where  [G]  is  now  a  composite  amplification  matrix,  [B]  and  [C] 
are  generic  implicit  and  explicit  matrices  respectively,  and 
the  constant  vector  of  boundary  conditions  and  zeros  has  been 
omitted.  It  is  apparent  then  that  the  amplification  matrix 
determines  the  numerical  solution  after  repeated  matrix 
multiplication  for  each  time  step  taken,  such  that: 

<U(i ,j ,n)>  =  lGJn<U(i, j,0)> 

A  second  way  to  obtain  the  amplification  matrix  describing  the 
FDM  is  to  use  the  Pade ’  Approach. 

Pade ’  Approach  to  the  Amplification  Matrix 

The  Pade’  approach  to  deriving  the  amplification  matrix 
is  presented  in  detail  by  Smith  (13:111-119).  This  approach 
begins  by  replacing  the  partial  differential  equation  with  a 
first  order  system  of  ordinary  differential  equations,  which 
are  then  solved  to  get  an  exponential  matrix  solution.  The 
exponential  matrix  solution  can  then  be  written  in  the  form  of 


a  recurrence  relation. 


Finally,  the  exponential  matrix 


<U ( t+k) >  =  exp(k [A] ) <U( t ) > 


or,  if  the  solution  vector  is  replaced  by  the  initial 
perturbation  vector  <E(0)>,  then  (13:113): 


<E( t+k) >  =  exp(k [A] ) <E(t ) > 


The  Pade ’  approximants  to  an  exponential  function  with  a 
real  argument  are  then  used  to  derive  a  rational  function, 
Rs.t(x),  which  is  called  the  (S,T)  Pade’  approximant ,  and  has 
the  form  (13:116): 


Rs.t(x)  =  (l+pix+. . .+pTxT )/(i+qlX+. . .+qsxs ) 

=  (PT (x ) } /{Qs (x) } 


Listed  below  are  three  Pade’  approximants  to  the  exponential 
matrix  that  result  in  three  well-known  FDMs  when  x  is  replaced 
by  k [ A]  (13:117). 


(S,T) 

Rs . t (x ) 

FDM 

(0,1) 

1+x 

Fully 

Expl icit 

(1.0) 

l/(l-x) 

Fully 

Impl icit 

(1.1) 

(l+x/2)/( l-x/2 ) 

Crank- 

■Nicol son 

The  recurrence  relation  for  the  EX  method  then  becomes 
equal  to  <E(t+k)>  =  ( [ I ] +k [A] ) <E( t ) >  ,  where  [I]  is  the 

identity  matrix.  So  for  n  time  steps  from  the  initial  error 
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vector,  <E(tn)>  =  (  [  I  ] +k  [  A]  )*»  < E  ( 0  )  >  .  The  recurrence  relation 
for  the  one  dimensional  case  is  also  the  recurrence  relation 
for  the  two  dimension  case  if  [A]  =  IB]  +  [C]  ,  where  [B] 
and  [C]  are  the  component  matrices  of  a  two  dimensional  FDM 
matrix  broken  up  into  block  matrix  form  (14:221).  Finally, 
the  amplification  matrix,  (G] ,  is  then  obtained  directly  from 
the  (S,T)  Pade ’  approximant  (13:122). 

Once  the  amplification  matrix  [G]  has  been  obtained,  it 
becomes  necessary  to  derive  all  the  eigenvalues  of  the 
amplification  matrix.  Smith  gives  the  general  equation  for 
the  eigenvalues,  Zs ,  of  a  common  tridiagonal  matrix  as 

, -  S7l 

Zg=b  +  2Yac  coS|q-  ,  s  =  1  (1)N  —  1  (13) 


where  a,  b,  and  c  are  the  three  elements  of  a  tridiagonal 
matrix,  from  left  to  right  (13:59).  For  the  standard  Pade' 
tridiagonal  matrix  elements,  the  eigenvalues  of  matrix  [A] 
then  become: 


z, 


4  .  2  S7t 


s  =  1  ( 1 )  N  —  1 


(14) 


For  two  dimensions,  the  eigenvalues  of  matrix  IA]  can  be 
derived  by  using  Equation  (13)  and  the  definition  of  (A)  =  (B] 


+  (C]  as  (14:221): 


Zmji= — 2\sin  2N+sin  2N  I  .  m=l(l)N-l,  n  =  1  (1)  N  -  1  (15 

The  eigenvalues  for  [A]  are  then  plugged  into  the 
equation  of  the  amplification  matrix,  (GJ ,  replacing  the 
matrix  [A],  as  shown  by  Smith  (13:120-124)  and  Twizell 
(14:212-215),  to  obtain  the  eigenvalues  of  the  amplification 
matrix,  (G] .  Substitution  of  the  N-l  unique  eigenvalues  is 
possible  because  the  associated  N-l  linearly  independent 
eigenvectors  can  be  used  as  a  basis,  or  coordinate  system,  in 
N-l  dimensional  space.  The  vector  of  initial  errors,  <E(0)>, 
can  then  be  represented  as  a  linear  combination  of  constants, 
Cs ,  multiplied  by  the  N-l  respective  eigenvectors  of  matrix 
(A).  Once  this  is  assumed  to  have  been  done,  the  eigenvalues 
operating  on  the  eigenvector  components  of  the  initial  error 
vector  replaces  the  matrix  multiplication  of  [A]<E(0)>.  This 
replacement  of  the  amplification  matrix  by  its  eigenvalues  is 
possible  because  of  the  definition  of  a  matrix  eigenvalue  and 
eigenvector  as  satisfying  (13:120): 

(A) <VS >=ZS <VS >  (16 

where  Zs  is  an  eigenvalue,  and  Vs  is  its  associated 
eigenvector.  Therefore,  the  error  grows  as  (13:120): 


<E(tn )> 


cs  [Rs.T(kZs)J“<Vs> 


s  =  l ( 1  )N-1  (17 


where  the  error  will  decay  only  if  |Rs.t(RZs)|  <  1  (13:120) 

Note  that  Rs,i(kZs)  represents  the  amplification  matrix 
eigenvalues,  and  may  be  referred  to  as  Gs ,  the  growth  factor 
associated  with  <Vs>.  The  amplification  matrix  eigenvalues, 
Gs  and  Gm ,  n ,  for  the  finite  difference  methods  used  in  this 
study  are  listed  below: 


EX1-D:  G,  -  1+kZ, 


EX  2-0:  Gmn=  l+kZ„ 


DF1-D:  G. 


2r  C,  +/-  V  1  -4r  (l  -C,j 


DF2-D:  Gm»  = 


-  * CC”  + ' c  J */- V  *  -  4r  l4  ~  (A  *  C.J2 


1  +0.5  kZ, 
CN  1-0:  G,  =  i  _o.5  kZ, 


CN  2-D:  Gmji = 


1  +0.5  kZ„ 


mjl  1  -  0.5  kZ_ 


ADI  2-D:  Gnvn  = 


r\  +0.5  kZ„Vl  +0.5  i 


(1-0.5  kZm)  (1-0.5  kZn) 


IM  1-D:  G,  = 


1-kZ. 


IM  2-D;  Gm,n  |  __  yy 


A  V  *  V vyiAA, A  *A  A  A  A  -rV,.  A  /.  A  .A  A  Ar  A  A  A  A  A  A  A  %A  A  - 


r»,»wn 


c 


where  Zs  =  Z*  =  Zn  =  Equation  (14),  Zm , n  =  Equation  (15), 

=  CB  =  Cn  =  cos  ( s  7T/N)  ,  and  the  indices  s>  ■ >  and  n  all 

vary  independently  as  1(1)N-1.  The  1-D  growth  factors  listed 
above  can  be  found  in  Smith  (13:120-123,150-53). 

Amplification  matrix  eigenvalues  for  the  2-D  cases  were 
derived  for  this  thesis  in  a  manner  similar  to  the  1-D  cases 
(see  Appendix  C  for  example). 

By  applying  the  same  criteria  to  the  maximum  and  minimum 
amplification  matrix  eigenvalues  as  was  applied  in  the 
coefficient  method  toward  the  p-ratios,  the  same  set  of  four 
FDM  stability  and  oscillatory  behaviors  can  be  predicted. 

Using  approximations  for  large  N,  several  predictions  based  on 
the  full  matrix  method  have  been  summarized  in  Figure  1.  It 
should  be  emphasized  that  the  matrix  stability  predictions  in 
Figure  1  are  the  full  matrix  method  predictions  based  on 
extreme  eigenvalues  and  large  N  approximations,  and  are 
therefore  worst  case  predictions  that  should  never  be 
violated.  Consequently,  it  is  not  surprising  that  the  full 
matrix  method  is  found  to  be  much  more  severe  in  its 
constraints  on  the  values  of  r,  requiring  smaller  values  of  r, 
than  the  coefficient  method. 

Non-Pi sastrous  Consequences  from  Oscillatory  Terms 

In  addition  to  predicting  the  stability  and  oscillatory 
behavior  of  FDMs ,  the  matrix  method  can  even  be  applied  to 
determine  when  the  oscillations  are  disastrous  or  not 
disastrous  to  the  solution.  This  is  extremely  useful  since 
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the  matrix  method  imposes  severe  constraints,  more  severe  than 
usually  needed,  because  the  range  of  eigenvalues  of  the 
amplification  matrix  vary  and  are  not  all  as  large  as  the 
extreme  values  used  in  the  predictions.  Lawson  and  Morris 
(10:1214)  pointed  out  that  oscillatory  terms  will  not  have 
disastrous  effects  if  the  highest-f requency  component, 
Cn-i<Vn-i>  ,  decays  to  zero  faster  than  the  lowest-f requency 
component,  Ci <Vi >  ,  where  from  the  equation  of  the  initial 
errors : 


<E ( 0 )  > 


N-l 

I 


s*l 


CS<VS> 


s= 1 ( 1 ) N-l 


in  which  E(0)  was  represented  as  a  linear  combination  of 
constants,  Cs,  multiplied  by  the  eigenvectors  of  the  [A] 
matrix  (10:1214).  Functional  relationships  for  the  above 
inequalities  have  been  derived  for  the  EX,  CN ,  and  IM  methods, 
and  are  listed  below  for  the  EX  and  CN  methods: 


1-D 

EX: 

2/k  - 

4/h2  >  7 T2/L2 

(18) 

2-D 

EX: 

2/k  - 

8/h2  >  27T2/L2 

(19) 

1-D 

CN: 

k/h  < 

L/7T 

(20) 

2-D 

CN: 

k/h  < 

l/(27T) 

(21) 

where  k  =  At,  h  =  Ax  »  Ay,  and  L  is  the  length  of  the  rod  or 
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the  length  of  one  side  of  the  square  FDM  mesh,  such  that 
N=L/h.  Also,  the  above  equations  are  approximations  based  on 
large  N. 

By  transforming  from  an  h,  k,  L  coordinate  system  to  an 
r,  N  coordinate  system  (see  Figure  2),  the  relationships  are 


easier  to  see  in  terms  of  the  familiar  r  coefficient.  The 
transformed  equations  are  listed  below: 


1- D  EX:  r  <  1 / { 2+7T 2 / ( 2N2 ) } 

2- D  EX :  r  <  l/(4+7T2 /N2 ) 

1- D  CN:  r  <  N/7T 

2- D  CN:  r  <  N/(27T) 


The  above  functional  relations,  where  r  is  set  equal  to 


the  expression  on  the  right  side  of  the  equations,  are  plotted 
on  the  h-k  domain  in  Figure  2.  These  restrictions  on  the 
functional  relationship  in  the  h-k  domain  are  less  severe  than 
the  stable  decay  restrictions  based  on  the  extreme  eigenvalues 
of  the  matrix  method. 


Relating  Coefficient  and  Matrix  Methods 

In  deriving  the  full  matrix  method  stability  predictions 
summarized  in  Figure  1,  extreme  eigenvalues  for  large  N 
approximations  were  used.  It  would  be  simpler  to  apply  a 
single  value  representing  all  the  eigenvalues  of  the  matrix 
method,  that  didn’t  depend  on  large  N  approximations,  in 
deriving  the  matrix  method  predictions.  Further,  in  picking  a 
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w 


single  value,  it  should  provide  good  predictions  if  the  single 
quantity  expresses  the  statistical  expectation,  or  arithmetic 
mean,  value  of  the  amplification  matrix  eigenvalues.  These 
expectation  values  are  discrete,  and  are  defined  for  Equation 
(14),  a  common  component  of  the  amplification  matrix 
eigenvalue  equations,  as  follows: 

1  N_1  2 

<Z*>  =  NZT  X  Z.  =  “~2  •  N^2  s  =  1,  2, ...  N  -  1 

*=  l  h 

where  N  is  the  number  of  spatial  node  points,  and  h  is  the 
spatial  step  size.  Using  this  result,  the  expectation  value 
of  the  1-D  EX  amplification  matrix  eigenvalues  is: 

EX  1-D:  <GS>  *  l-2r 

where  for  this  equation,  <Gs>  denotes  the  expectation  value 
of  G$  . 

Comparing  the  1-D  EX  result  above  with  the  1-D  EX 
coefficient  method  p-ratio,  an  interesting  relationship  is 
seen  .  .  .  both  equations  are  the  same.  In  effect,  the 
coefficient  method  has  been  derived  from  the  more 
mathematically  rigorous  matrix  method.  This  is  an  important 
result  since  it  indicates  that  the  coefficient  method  is  not 
simply  a  crude,  single  node  point  approximation  applied  to  all 
the  node  points.  Instead,  the  coefficient  method  is  firmly 


related  to  the  matrix  method  and  statistically  considers  all 


the  nodes. 

Next,  the  above  results  and  conclusions  need  to  be 
checked  with  the  other  FDMs  used  in  this  study.  For  the  FDMs 
in  this  study,  the  expectation  values  for  all  the  1-D  and  2-D 
amplification  matrix  eigenvalues  were  derived  and  are 
tabulated  in  Appendix  D.  The  results  of  the  derivations 
indicate  that  the  expectation  value  of  the  amplification 
matrix  eigenvalues  is  equal  to  the  coefficient  method  p-ratio 
for  two-level  FDMs.  For  three-level  FDMs,  the  expectation 
values  diverge  from  the  p-ratios,  as  would  be  expected,  since 
the  behavior  of  the  three  FDMs  in  this  study  which  have  the 
same  p-ratio  is  obviously  different.  Therefore,  the 
coefficient  method  is  seen  to  be  a  worse  approximation  for 
three-level  FDMs,  but  may  still  be  statistically  related  to 


The  computational  approach  used  to  examine  stability 


predictions  will  be  discussed  in  this  chapter.  Topics 
included  will  be  the  domain  approach,  stability  components, 
how  to  read  the  domain  graphs,  and  comparisons  of  domain 
graphs  with  representative  conventional  graphs. 

Domain  Graphs 

Domain  graphs  provide  a  simple  and  comprehensive  way  to 
compare  all  functional  relationships  based  on  h,  k,  L,  and 
time  for  predicting  FDM  stability,  oscillations,  and  accuracy. 
Because  it  is  the  most  basic  of  all  possible  representations, 
the  h-k  domain  graph  is  one  of  the  best  ways  to  display  a  vast 
array  of  possible  coefficients.  Almost  all  other  coefficients 
and  conventional  graphs  can  be  found  as  simple  points  or  lines 
on  a  domain  graph.  Another  advantage  of  domain  graphs  is  that 
they  are  extendible.  Therefore,  the  predicted  functional 
relationship  does  not  have  to  be  known  and  built  into  the 
graph  when  it  is  made.  Consequently,  the  domain  graph  is 
useful  for  comparing  other  stability  predictions  in  the 
future.  However,  the  domain  graph  does  have  one  chief 
drawback.  Limited  domains  are  a  problem  since  the  number  of 
computations  needed  to  solve  many  heat  equations  increases 
approximately  as  the  total  number  of  node  points  that  must  be 
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calculated  in  the  domain.  For  large  numbers  of  maximum  time 
and  space  steps  on  the  domain,  the  total  number  of  node  points 
to  be  calculated  approaches: 


1-D, 

genera  1 : 

(nx2/2) (nt 2/2) 

and 

(22) 

2-D , 

Ax  =Ay  plane: 

(nxv3/3) (nt  2/2) 

and 

(23) 

2-D, 

genera  1 : 

(nx2/2)(ny2/2)(nt2 

/  2 ) 

(24) 

where  nx 

and  riv  are  the 

maximum  number  of 

s  t  eps 

i  n 

the  x  and  y 

directions  respectively 

,  nx y  is  the  maximum  number 

of  steps  in 

the  x  and 

y  directions 

for  the  Ax=Ay  plane 

(a  square  grid  is 

assumed , 

so  x  and  y  steps  are  equal),  and 

nt  is 

the 

maximum 

number  of 

steps  in  time 

to  reach  the  solution. 

In 

addit ion , 

all  steps  up  to  and  including  the  above  maximum  numbers  are 
assumed  to  be  used,  starting  with  one  step  in  each  dimension. 
To  simplify  comparison  of  the  domain  graphs,  the  same  domains 
were  used  on  all  graphs.  As  a  result,  equations  (22)  and  (23) 
show  that  the  domain  graphs  in  this  study  were  limited  by  the 
two  dimensional  case.  Ultimately,  the  domain  graph  is  a 
technique  that,  with  the  advent  of  faster  and  more  powerful 
digital  computers,  promises  to  see  more  use  in  the  analysis  of 
methods  for  predicting  FDM  behavior. 

Sta bi 1 i ty  and  Behav i ora  1  Components 

After  having  selected  the  h-k  domain  as  the  best  way  to 
compare  various  stability  properties,  oscillations,  and  final 
errors,  one  of  the  major  difficulties  that  had  to  be  overcome 


in  using  the  approach  was  that  of  precisely  defining  what  was 
to  be  displayed  in  the  h-k  domain,  and  then  designing 
algorithms  that  allowed  the  computer  to  automatically  find  and 
quantify  the  information  to  be  displayed.  Various  FDM 
behavioral  diagnostics  were  examined,  and  the  following  were 
selected  as  a  minimal,  yet  complete,  set:  maximum  number  of 
time  oscillations,  maximum  number  of  spatiai  error 
oscillations,  growth  or  decay  of  errors,  and  maximum  error  in 
the  final  solution. 

In  chapter  3,  stability  analysis  using  the  p-ratio  and 
the  eigenvalues  of  the  matrix  method  predicted  four  FDM 
behaviors.  These  four  behaviors  can  concisely  be  displayed  as 
combinations  of  individual  stability  components  separated  onto 
two  graphs:  one  graph  for  the  growth  or  decay  of  roundoff 
erroi  as  determined  by  the  whether  the  numerical  solution  is 
bounded,  and  a  second  graph  for  showing  whether  the  numerical 
solution  oscillates  or  not. 

In  a  previous  paper  (8:29),  the  basic  method  used  to 
prove  stability  was  to  show  decreasing  rms  and  maximum  error 
as  functions  of  time  for  a  given  solution  of  the  heat 
equation.  This  diagnostic  was  adopted,  but  since  stability  is 
only  a  necessary  condition  for  accuracy,  this  diagnostic  can 
only  be  used  to  prove  stability,  and  can  not  be  used  to  prove 
a  solution  is  unstable.  So  if  the  numeric  versus  analytic 
solution  errors  grow  in  time,  then  the  solution  may  or  may  not 
be  stable  with  respect  to  the  numerical  roundoff  errors. 
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However,  since  all  the  methods  are  unconditionally  stable, 
except  for  the  EX  method,  these  graphs  play  a  minor  role  in 
this  thesis,  and  have  only  been  included  in  Appendix  A  to  show 
that  they  are  an  unreliable  stability  diagnostic.  To 
implement  this  diagnostic  so  that  a  digital  computer  could 
automatically  determine  the  trend  of  error  growth  with  time,  a 
line  was  curve  fitted,  using  least  squares  techniques,  to  the 
maximum  error  versus  time  at  each  time  step.  Finally,  if  the 
slope  of  the  fitted  line  is  negative,  indicating  decreasing 
error,  the  FDM  is  stable  at  that  point  on  the  h-k  domain. 

Since  stability  is  a  concept  based  on  the  growth  or  decay 
of  numerically  induced  errors  and  oscillations  in  time,  all 
nodes  for  each  heat  equation  solution  in  the  h-k  domain  were 
tracked  with  respect  to  time.  The  maximum  number  of  slope 
changes  for  any  of  the  single,  fixed  spatial  nodes  was  counted 
as  the  number  of  oscillations  for  stability  purposes.  In 
addition,  spatial  oscillations,  although  not  directly 
predicted  by  p-ratios  or  matrix  eigenvalues,  are  graphed.  Two 
reasons  for  including  spatial  oscillations  exist.  One, 
various  texts  (13:124;  14:206-207,213)  show  spatial 

oscillations  when  discussing  stability,  and  the  danger  could 
occur  of  associating  stability  predictions  with  spatial 
oscillations  in  the  final  solution,  as  was  incorrectly  done  in 
a  previous  paper  (8:29).  Secondly,  oscillations  predicted  by 
stability  methods  in  time  are,  for  practical  purposes,  only 
really  important  if  they  result  in  spatial  oscillations  in  the 
final  solution.  Since  the  actual  solution  in  a  general 
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partial  differential  equation  could  be  oscillating,  the 
analytic  solution  was  subtracted  from  the  numerical  FDM 
solution  to  get  the  actual  number  of  slope  change  oscillations 
with  respect  to  the  analytic  solution.  So,  the  spatial 
oscillation  graphs  display  the  maximum  number  of  error  slope 
change  oscillations  in  either  spatial  direction  for  the  two 
dimensional  FDMs . 

Lastly,  whether  a  solution  oscillates  or  not  is  a 
secondary  consideration  to  the  final  maximum  error  in  the 
solution.  Also,  since  the  matrix  method  can  even  be  used  to 
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generate  functional  relationships  in  the  h-k  domain  that 
predict  when  the  oscillations  are  disastrous  or  not  disastrous 
to  the  final  solution,  these  predictions  must  be  examined. 
Therefore,  h-k  domain  graphs  of  the  maximum  error  of  all  the 
node  points  in  the  final  solution  have  been  included  in  the 
analysis.  For  all  error  calculations,  the  error  used  is  the 
absolute  error,  defined  as:  error  =  1 U ( i , j , n ) -T( x , y , t ) I  , 
where  U  is  the  numerical  solution  to  the  FDM  and  T  is  the 
analytic  solution  to  the  partial  differential  equation. 

Reading  h-k  Domain  Graphs 

The  domain  graphs  used  in  this  thesis  are  not 
conventional,  so  this  section  provides  the  extra  information 
needed  to  make  using  the  domain  graphs  easier  and  more 
effective.  First,  two  partial  differential  equation  problems 
were  used,  one  for  the  one  dimensional  case,  and  one  for  the 
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two  dimensional  case.  The  domains  were  all  chosen  to  cover 


the  same  range  in  spatial  and  temporal  step  sizes,  so  that  the 
resulting  graphs  could  be  compared  directly.  Logarithmic 
scales  were  used  on  the  axes  because  most  of  the  predicted 
functional  relationships  become  simple  lines  on  a  logarithmic 
scale.  Also,  since  exact  solutions  are  computed,  the  h-k 
domain  is  made  up  of  a  discrete  set  of  step  sizes  that  exactly 
step  through  the  space  and  time  dimensions  in  nx ,  ny ,  and  nt 
steps.  Therefore,  the  values  of  the  z-axis  are  actually  only 
computed  at  discrete  points,  even  though  they  are  shown  as  a 
continuum  of  values.  As  a  result,  the  contour  graphs 
introduce  winding  ribbons  and  waves  that  are  mainly  artifacts 
of  the  discrete  points  used.  Reference  bottom  Figure  3,  where 
lines  intersect,  to  see  the  locations  of  the  2-D  discrete 
points. 

The  graphs  for  both  the  one  and  two  dimensional  transient 
heat  equations  FDMs  are  all  three  dimensional  graphs  plotted 
as  two  dimensional  contour  graphs,  and  are  read  very  much  like 
two  dimensional  contour  maps  used  in  navigation.  These  graphs 
can  be  broken  down  into  just  four  types,  depending  on  what  the 
z-axis,  which  is  represented  by  the  contour  lines,  itself 
represents.  These  graphs  were  chosen,  because  while  harder  to 
visualize  topography  on  than  three  dimensional  graphs,  they 
are  much  easier  to  use  to  visualize  the  predicted  functional 
relationships  on  than  three  dimensional  graphics.  The  contour 
lines  can  then  represent  the  following:  1)  maximum  numeric 

node  time  oscillations,  2)  maximum  spatial  error  oscillations, 


m* w. 


3)  maximum  absolute  error  in  the  final  solution,  and  4)  the 
slope  of  the  maximum  absolute  error  versus  time  (found  in 
Appendix  A).  Otherwise,  all  horizontal  and  vertical  axes  are 
the  same  for  all  the  graphs,  all  the  h-k  domains  are  the  same 
for  each  graph,  and  all  the  graphs  are  the  same  size  for  ease 
of  use  and  behavioral  comparison.  In  addition,  both  the  one 
and  two  dimensional  solutions  for  each  FDM  will  generally  be 
placed  together  on  a  page,  also  for  easy  comparison  of 
behavior  changes.  All  the  horizontal  axes  represent  the 
logarithms  of  the  x-steps.  In  the  two  dimensional  case,  the 
x-step  =  y-step  plane  is  used  for  3  reasons:  1)  most  FDMs  are 
normally  derived  for  this  plane,  2)  the  same  form  of  h-k 
domain  graph  can  be  used,  and  3)  it  reduces  the  total  number 
of  nodal  points  to  be  calculated  per  domain.  All  vertical 
axes  represent  the  logarithms  of  the  time  steps,  and  are 
exactly  the  same  for  both  the  one  and  two  dimensional  graphs. 

To  actually  read  a  z-axis  iso-value,  the  following  should 
be  observed.  Contour  line  numbers  increase  with  increasing 
values,  after  any  z-axis  mapping  transforms  have  been  used,  of 
the  z-axis.  Contour  line  numbers  are  always  located  to  the 
right  of  their  respective  contour  line.  The  actual  value  of 
the  z-axis  iso-level  line  can  be  obtained  by  using  the 
following  equation: 


( L#-SL# ) (Of f set )+( Start  z  value) 


where  L#  is  the  number  of  the  line  you  wish  to  find  the  z-axis 
value  of,  SL#  is  the  number  of  the  line  for  the  starting  z 
value,  and  Offset  is  the  +  /-  offset  value  used  for  each 
contour.  Once  equation  (25)  has  been  calculated,  the 
appropriate  inverse  z-axis  transform  must  be  performed.  Four 
z-axis  transforms  have  been  used:  linear,  loglO,  modloglO, 

and  mod21ogl0,  where: 

modloglO(z)  =  sgn ( z ) *  1 ogl 0 ( abs ( z ) +  1 ) 

mod2 1 ogl 0 ( z )  =  mod  1 oglO ( 1 ogl 0 ( z ) ) 

with  mod21ogl0  only  being  used  for  the  EX  method  to  display 
the  maximum  errors  in  the  final  solution,  which  are  all 
positive  z  values  of  extreme  dynamic  range. 

However,  it  is  not  the  intent  of  this  paper  to  require 
careful  reading  and  extraction  of  values  off  the  contour 
plots,  and  they  can  be  used  in  a  much  simpler  fashion.  To 
simplify  reading  important  information  off  the  contour  graphs, 
two  types  of  contour  lines  are  used,  solid  and  dotted.  For 
oscillation  graphs,  regions  of  the  h-k  domain  with  solid  lines 
are  oscillating  and  regions  of  the  domain  with  dotted  lines 
are  not  oscillating.  For  maximum  error  graphs,  dotted  contour 
lines  show  regions  where  the  maximum  solution  error  is  less 
than  10%  of  the  range  of  the  actual  analytic  solution’s 
values.  Therefore,  regions  of  a  maximum  error  graph’s  h-k 
domain  with  dotted  lines  have  maximum  errors  that  are  not  too 
disastrous  to  the  final  solution.  Finally,  for  the  slope  of 


the  maximum  error  graphs  in  Appendix  A,  negative  slopes 
indicate  stable  regions  and  are  shown  as  dotted  lines,  while 
the  stability  of  solid  lines  can  not  be  proved  using  this 
diagnostic . 

Comparison  of  h-k  Domain  Graphs  with  Conventional  Graphs 

Comparison  of  the  h-k  Domain  Graphs  with  conventional 
graphs  is  accomplished  for  two  reasons.  One,  it  provides 
verification  of  the  computer  implementations  of  the  FDMs  and 
diagnostics  using  previous  published  results.  Secondly,  it 
shows  how  to  understand,  conventionally,  what  the  h-k  domain 
graphs  represent.  Figures  4  and  5  show  the  one  dimensional 
analogs  for  points  on  the  h-k  domain  versus  conventional 
graphs,  where  the  number  of  time  steps  is  10,  and  space  steps 
is  14  (using  Problem  1  of  Chapter  5).  For  this  point  in  the 
h-k  domain,  the  diagnostics  predict  the  following  values: 


FDM 

time 

osc 

space 

osc 

max 

error 

slope 
max  error 

time 

r 

1-D 

IM: 

0 

1 

.  0329 

-0.206 

.25s 

4 . 90 

1-D 

CN  : 

8 

9 

.0104 

-1 . 09 

.  25s 

4 . 90 

where  the  underlined  values  can  be  found  on  Figures  4  and  5. 

The  1-D  CN  method  is  compared  side  by  side  with  the  1-D 
IM  method.  The  IM  method  is  used  as  a  reference  since  the  IM 
method  is  known  to  be  non-osc i 1 1 atory  due  to  its  Lo-stability 
As  seen  in  Figures  4  and  5,  the  program  correctly  identifies 
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the  maximum  error,  the  negative  slope,  and  the  number  of  slope 
changes  related  to  oscillation::.  Note  that  the  number  of 
spatial  oscillations  for  the  CN  case  is  nine,  and  not  the 
seven  easily  counted  on  the  graph,  since  the  rise  in  the 
center  of  the  graph  is  too  small  to  be  seen. 

Figures  6  and  7  use  graphs  from  two  previous  papers  on 


stabi 1 ity 

and 

the 

programs 

in  this  thesis 

have  been 

used  to 

calculate 

behavior 

that  should  be 

seen  in 

the  graphs 

(4:23,24; 

8:42, 46) . 

The  results  are 

as  foil ows : 

FDM 

time  space 
osc  osc 

max 

error 

slope  time 

max  error 

r 

2-D 

IM: 

0 

1 

8.06 

-298 

.  4s 

10.0 

2-D 

ADI  : 

2 

3 

10.5 

-1058 

.  4s 

10.0 

2-D 

CN: 

2 

5 

148 

-425 

.  4s 

10.0 

2-D 

CN: 

6 

7 

4.45 

-674 

.  2s 

2.5 

2-D 

DF : 

6 

1 

18.7 

-.0707 

.  2s 

1 . 0 

where  again, 

the  implicit 

method 

provides 

a  reference  value, 

and  the  underlined 

values 

can  be 

found  in 

the  figures.  Note 

that  the 

1  ast 

two 

sets  of 

numbers 

are  at 

,  2s  so  that 

the 

values  can  also  be  read  off  the  appropriate  domain  graphs  in 
this  thesis.  These  values  agree  well  with  the  previous  thesis 
graphs.  However,  the  previous  thesis  graph  in  Figure  7  only 
shows  the  oscillatory  behavior  in  time  of  one  spatial  node. 
Since  the  point  calculated  in  the  h-k  domain  is  the  maximum 
value  of  all  the  spatial  nodes,  the  one  node  graph  only  shows 


5  slope  change  oscillations,  and  not  the  maximum  value  of  six. 
Note  that  in  Figure  6,  the  CN  graph  of  the  maximum  error 
versus  time  for  a  fixed  solution  time  and  various  values  of  r 
would  trace  out  a  straight  line  on  the  h-k  domain,  and  shows 
the  power  and  versatility  of  being  able  to  represent  several 
types  of  line  graphs  for  several  types  cf  prediction 
coefficients  on  one  h-k  domain  graph.  Finally,  note  that  zero 
time  oscillations  result  from  the  implicit  method,  so  zero 
oscillations  denotes  no  time  oscillations,  yet  1  spatial  error 
oscillation  means  the  solution  is  not  oscillating  in  space. 
This  apparent  contradiction  occurs  because  the  solution  of  the 
analytic  solution  has  one  slope  change,  and  a  monoton i ca 1 1 y 
decreasing  IM  solution  also  has  one  slope  change.  Therefore, 
one  slope  change  is  not  considered  a  spatial  error  oscillation 
for  the  analytic  problems  used  in  this  thesis. 


c 
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V .  Results 

This  chapter  presents  the  h-k  domain  graphs  generated  by 
using  the  EX,  DF,  CN,  ADI,  and  IM  methods.  These  graphs  were 
generated  using  two  sample  problems  with  known  analytical 
solutions.  First,  the  two  problems  will  be  examined.  Then, 
the  h-k  domain  graphs  will  be  investigated  and  the  predictions 
from  the  coefficient  and  matrix  methods  compared  with  the 
graphs . 

Problem  #1 

The  first  problem  is  used  in  conjunction  with  all  the  one 
dimensional  FDMs ,  and  was  used  by  Smith  (13:124)  in  discussing 
stability,  and  is  equivalent  to  problems  used  by  Lawson  and 
Morris  (10:1212)  and  Twizell  (14:205-206)  in  studying 
stability.  The  problem  is: 

U ( x , t ) t  =  U ( x , t ) x x  0  <  x  <  1  t  >  0  satisfying 

U ( 0 , t )  =  U  ( 1  ,  t )  =  0  t>0 

U ( 0 , x )  =1  0  <  =  x  <  =  1  t  =  0  and 

so  that  x  varies  from  xi=0  to  X2=l  in  length,  and  time  t 
varies  from  ti=0  to  t2=0.25  in  seconds.  Partial 
differentiation  of  U(x,t)  with  respect  to  t  and  x  is  denoted, 
in  standard  concise  notation,  by  subscripted  variables.  For 
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the  1-D  CN  FDM,  tj  is  also  set  equal  to  0.50  seconds  to 
compare  the  domain  graphs  at  two  different  final  times  so  that 
the  domains  remain  constant,  but  the  number  of  time  steps  used 
increases.  The  general  analytic  solution  is  found  in  Carslaw 
and  Jaeger,  and  for  problem  #1  becomes  (2:96): 


U(x,t) 


E4  2  2 

—  exp(-n  7t  t)  sin  (nrcx) 


n  =  1  ,  3 ,  5 , 


Problem  #2 


The  second  problem  was  used  in  two  previous  theses  to 
discover  oscillatory  behavior  and  examine  the  accuracy  of  the 
EX,  DF,  CN,  ADI,  and  IM  FDMs  (A: 19, 8: 19).  The  problem  is: 


U ( x , y , t ) t  =  U(x,y,t)*x  +  U(x,y,t)yy 


with 


U ( x , 0 , t )  =  U ( x , 1 , t )  =  U ( 1 , y , t )  =  U ( 0 , y , t )  =  400 


U(x ,y , 0) 


where 


0  <  x  <  1 


and  0  <  y  <  1 


The  analytic  solution  is  (0:19-20) 


U ( x , y , t )  =  400  (l-f(x,y,t)) 


where 


f (x.y.t)  =  g(x,t)  h(y , t ) 


where 


i 
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g(x,t)  =  ^  exp(-n2n2t)  sin  (nnx) 


h(y,t)  *  ^  ^  ~  exp(-n2u2t)  sin  (nny) 


n  =  1 ,  3 ,  5  ... 


h-k  Domain  Graph  Solutions 


Approximately  25  and  1.6  million  total  node  points  were 


calculated  in  the  2-D  and  1-D  problem  respectively,  in  order 


to  calculate  each  set  of  FDM  domain  graphs.  The  2-D  problem, 


as  was  expected  from  Equations  (22)  and  (23),  limited  the  size 


of  the  domains  that  could  be  calculated,  since  the  VAX 


computer  could  calculate  an  average  of  1.1  million  node  points 


per  hour.  Double  precision  accuracy  was  maintained  with 


minimum  series  term  calculations  in  the  analytic  solutions. 


The  minimum  number  of  series  terms  needed  at  each  time  step 


was  calculated  in  the  programs  based  on  curve  fitted  power 


functions  approximating  the  minimum  terms  needed  versus 


solution  time.  Oscillations  are  given  as  the  number  of  slope 


changes  and  maximum  absolute  errors  are  given  in  degrees.  The 


slope  of  the  maximum  error  versus  time  has  dimensions  of 


degrees/second  (Appendix  A).  The  spatial  step  size  associated 


with  the  horizontal  axes  have  units  of  length,  while  the 


vertical  axes  representing  the  time  step  duration  have  units 


of  seconds.  The  specific  units  for  degrees  and  length, 
whether  °K  or  °F,  meters  or  feet,  will  depend  on  the  choice  of 
thermal  diffusivity  (a),  which  has  a  value  of  i  for  the 


problems  used  in  this  thesis. 

F  u 1 1 y _ Imp  licit  Results 

The  IM  method  was  included  in  this  study  as  a  reference 
for  comparing  the  behavior  of  other  FDMs ,  since  it  was  known 
to  be  Lo-stable  (14:214).  As  a  result,  no  IM  h-k  domain 
graphs  for  oscillations  have  been  included,  since  the  entire 
domains  are  devoid  of  contour  lines.  No  contour  lines  can  be 
drawn  because  all  the  points  in  the  h-k  domain  return  0  as  the 
maximum  number  of  numeric  node  time  oscillations,  and  1  as  the 
maximum  number  of  spatial  error  oscillations. 

The  IM  maximum  error  graphs  show  the  least  influence  of 
an  r  type  functional  dependence  in  the  domain.  From  Figure  8, 
the  accuracy  of  the  IM  method  appears  to  depend  largely  on 
using  smaller  time  steps,  with  less  dependence  on  the  space 
step  size,  once  a  certain  number  of  space  steps  has  been  used. 

Fully  Explicit  Results 

The  results  of  the  EX  method  provide  easily  understood 
graphs  in  the  h-k  domain,  and  are  presented  in  Figures  9 
through  11.  The  r  correlation  in  the  graphs  is  extremely 
high,  with  the  solid  contour  lines  representing  time 
oscillations  winding  back  and  forth  about  the  r=0.5  line  for 
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the  1-D  case,  and  r=0.25  for  the  2-D  case.  Note  that  the 
discrete  winding  back  and  forth  of  the  contour  lines  appears 
to  be  largely  an  artifact  of  using  discrete  solution  points 
(see  Figure  3).  The  best  fit  r  values  mentioned  above  match 
exactly  with  the  predicted  coefficient  method  values  where 
r>0.5  and  r>0.25  are  predicted  as  regions  where  time 
oscillations  begin  to  occur.  As  a  matter  of  fact,  it  is  seen 
from  Figure  9  that  as  r  increases  from  the  predicted  starting 
values  for  the  oscillations,  a  sort  of  turbulent  zone  is 
encountered  where  the  number  of  oscillations  increases  until 
it  reaches  the  maximum  number  of  slope  changes  that  the  number 
of  time  step  nodes  will  allow,  and  then  remains  at  that 
maximum,  as  shown  by  the  horizontal  lines  reaching  to  the  time 
axis.  Therefore,  a  line  perpendicular  to  the  r  lines  would 
show  that  with  increasing  r  the  oscillations  also  increase 
until  the  maximum  time  oscillations  allowed  are  reached.  The 
matrix  method  predictions  for  the  start  of  1-D  and  2-D 
oscillations  are  r>0.25  and  r>0.125  respectively.  So  the 
matrix  method  is  better  at  predicting  when  no  oscillations 
will  occur,  while  the  coefficient  method  appears  better  at 
predicting  when  oscillations  start  to  occur,  for  these  example 
problems.  The  actual  EX  method  r  variation  is  small  compared 
to  the  coefficient  predictions,  and  by  itself  is  not 
significant  enough  to  discredit  the  coefficient  method  as 
generally  valid. 

Figure  10  shows  the  corresponding  results  for  spatial 
error  oscillations.  In  this  case,  uncharacteri st i c  of  the 


rest  of  the  FDMs  examined,  the  spatial  oscillations  also  match 
the  same  r  values  as  the  time  oscillations  in  Figure  9,  so 


predicted  stability  oscillations  in  time  correlate  well  with 
the  actual  spatial  oscillations  in  the  final  solution.  Also, 
the  same  type  of  turbulence  is  seen,  but  this  time  the 
oscillations  increase  until  the  maximum  number  of  spatial 
oscillations  is  reached,  as  shown  by  the  vertical  lines. 

Maximum  errors  for  each  solution  point  are  graphed  in 
Figure  11.  Once  again,  these  match,  uncharacteristically, 
very  closely  with  the  r  values  observed  for  the  oscillations. 
Based  on  the  matrix  method  error  predictions  for  large  N,  the 
matrix  method  is  successful  at  predicting  what  r  lines  below 
which  the  errors  from  the  oscillations  are  not  disastrous, 
even  when  N  is  only  about  3  on  the  domain  graph.  This  is  very 
good  for  a  first  order  approximation  of  the  trigonometric 
functions  in  the  eigenvalues  of  the  amplification  matrix. 

PuFort-Franke 1  Results 

The  DF  results  involve  jump  starting  the  three-level 
method  with  three  different  methods  to  get  the  first  time 
level.  Both  the  analytic  solution  and  IM  jump  starts  are 
graphed  for  both  the  1-D  and  2-D  cases,  and  the  CN  method  is 
used  to  jump  start  the  1-D  method  as  well.  Figures  12-14  and 
15-17  are  the  results  of  the  analytic  solution  and  IM  jump 
starts,  respectively,  while  Figures  18-19  show  the  results  of 
using  the  CN  jump  start  to  the  DF  method. 
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Figure  14.  DF  Maximum  Errors,  Analytic  Start 
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Figure  17.  DF  Maximum  Errors,  IM  Start 


The  time  oscillatory  behavior  displayed  by  the  DF  method 
is  much  more  complex  than  that  of  the  EX  method.  However,  all 
the  time  oscillation  graphs  show  high  r  correlations  in  the 
slopes  of  the  lines  best  fitting  where  the  oscillations  begin. 
But,  unlike  the  EX  method,  the  number  of  oscillations  at  first 
increases  with  increasing  r  values,  and  then  begins  to 
decrease  again,  and  even  stops  oscillating  for  larger  values 
of  r.  This  is  a  remarkable  result,  and  should  be  looked  into 
in  more  detail.  Unfortunately,  the  full  matrix  method 
eigenvalues  are  complicated,  even  for  the  1-D  case,  since  for 
values  of  r  greater  than  about  0.5,  the  eigenvalues  generally 
become  complex  valued.  Due  to  the  complex  eigenvalues  of  the 
amplification  matrix  for  the  1-D  and  2-D  DF ,  a  complete 
analysis  of  the  DF  method  was  not  completed.  Furthermore, 
complications  can  arise  when  h  is  approximately  equal  to  k 
because  "the  computed  solution  may  involve  waves  associated 
with  the  hyperbolic  character  of  the  PDE  (4.4.3),  with  which 
it  is  consistent"  (9:167).  This  hyperbolic  wave  region  of  the 
h-k  domain  also  includes  the  region  observed  to  have  temporal 
and  spatial  oscillations  decreasing  to  zero.  Waves  in  the 
jump  start  solution  may  also  be  the  reason  that  the  CN  jump 
start  has  the  lowest  resulting  oscillatory  behavior.  The 
reason  for  this  behavior  needs  further  study,  since  the 
observed  behavior  would  be  expected  to  be  just  the  opposite, 
with  the  analytic  having  the  least,  the  IM  more,  and  the  CN 
jump  start  the  most  oscillations  for  the  same  h-k  domain. 
However,  since  the  errors  for  such  large  values  of  r  are 


generally  increasing,  even  though  no  oscillations  may  be 

#  occurring,  the  solutions  are  unusable.  In  fact,  the  solutions 
are  about  2  orders  of  magnitude  worse  that  the  IM  method, 
which  doesn’t  oscillate  in  high  r  regions  either. 

(9  For  the  1-D  and  2-D  DF ,  the  coefficient  method  predicted 

r>0.5  and  r>0.25  respectively,  for  the  start  of  oscillatory 
behavior.  The  domain  graphs  show  that  the  actual  values  of  r 

#  are  between  0.50  to  1.0  and  between  0.25  and  0.50  for  the  1-D 
and  2-D  results  respectively.  However,  the  analytic  1-D  jump 
start  takes  exception,  with  the  major  r  correlation  best  fit 

f*  line  for  the  start  of  oscillations  falling  at  slightly  less 

than  r=0.5.  Also,  for  the  1-D  CN  jump  start,  the  best  fit  r 
line  moves  to  about  r=1.0.  As  such,  the  coefficient  method 

#  predicted  behavior  is  generally  good  at  predicting  when  the 
oscillations  do  not  occur,  but  doesn’t  do  as  well  as 
predicting  where  the  oscillations  actually  begin.  This  is 

#  expected,  since  for  any  given  test  problems,  the  behavior  of 
the  FDN  should  be  better  and  not  worse  than  the  worst 
predicted  behavior.  Finally,  the  coefficient  method  is 
inadequate  as  a  tool  for  explaining  any  of  the  complex 
behavior  seen  in  the  DF ,  while  the  matrix  method  may  provide 
some  i ns i ght  s . 

The  spatial  oscillations  for  the  DF  generally  have  a  more 
linear  power  function  form,  since  they  have  less  slope  than 
the  r  coefficient  lines  designating  the  start  of  osci 1 lations . 
Two  other  results  are  worth  noting.  One,  the  1-D  case  for  the 
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analytic  and  IM  jump  starts  oscillate  spatially  throughout 
most  of  the  domain,  except  in  the  high  r  region  where  h  is 
approximately  equal  to  k.  Secondly,  the  1-D  DF  CN  jump  start 
and  both  of  the  2-D  DF  jump  starts  exhibit  a  more  linear 
functional  constraint  on  where  oscillations  start.  Plus,  the 
domain  graphs  look  similar,  with  low  peaks  of  oscillations 
occurring  at  locations  in  the  domain  containing  least  maximum 
error  potholes.  This  is  suspicious  of  a  lack  of  accuracy  and 
precision,  or  both,  creating  artificial  error  oscillations  in 
the  final  solution.  None  the  less,  it  seems  these  are 
actually  oscillations.  Because  the  smallest  maximum  errors  in 
the  domain  involve  5  to  7  digits  of  precision,  while  the 
program  should  be  accurate  to  at  least  10  digits  of  precision, 
neither  accuracy  nor  precision  seem  to  be  the  cause. 
Regardless,  the  oscillations  occur  at  minimum  potholes  of  the 
maximum  error,  and  can  be  completely  ignored  from  a  practical 
engineering  standpoint.  Concluding  spatial  oscillations,  the 
correlation  between  time  oscillations  predicted  by  stability 
and  actual  spatial  oscillations  is  shown  to  be  a  poor  one,  and 
that  spatial  oscillations  should  not  be  used  to  validate 
stability  oscillation  predictions,  as  is  often  inferred. 

The  maximum  error  domain  graphs  show  that  the  least 
errors  in  the  domain  generally  lie  in  valleys  with  high  r 
correlations,  where  water  would  flow  from  the  domain  if  it 
were  three  dimensional,  and  with  the  r  coefficient  becoming 
more  linear  as  the  10%  solid  contour  lines  are  reached.  The 
valleys  for  the  1-D  DF  are  approximately  centered  between 


r=l/8  to  1/2  (1/4  average)  for  the  1-D  cases,  and  between 
r=l/l6  to  1/4  (1/8  average)  for  the  2-D  cases.  The  analytic 
jump  starts  resulted  in  the  smallest,  most  severe,  r  values. 

Cra nk-Kicolson  Results 

Figures  21-23  show  the  standard  problems  at  the  standard 
times  for  the  one  and  two  dimensional  cases.  Figures  19-20 
also  show  the  1-D  CN  results,  but  at  a  final  solution  time  of 
0.50  seconds,  instead  of  the  standard  0.25  seconds.  This  was 
done  to  see  if  the  characteristics  of  the  domain  graphs  would 
change,  and  they  did  not.  The  similar  characteristics  for  the 
different  times  were  expected.  Since  the  domains  were  kept 
constant,  the  predictions  remained  constant  as  well.  The  only 
major  variation  was  that  the  time=0.50  second  domains  required 
using  twice  as  many  time  steps  to  create  the  same  domain 
graph.  Therefore  the  maximum  number  of  time  oscillations  that 
can  be  supported  by  the  nodes  also  doubled,  increasing  the 
turbulent  zone  area  required  to  reach  maximum  time 
oscillations  slightly,  but  still  starting  at  exactly  the  same 
r  value  centered  at  r=2.0.  Whereas  the  spatial  oscillation 
line  is  parallel  for  both  times,  it  become  less  constraining 
and  moved  up  higher  in  the  domain  for  the  t=0.50  second  case. 
The  domain  itself  supports  about  the  same  number  of  spatial 
oscillations,  though  a  slight  decrease  was  observed.  The  10% 
maximum  error  line  falls  exactly  in  the  same  place  for  both 
cases.  However,  the  valley  (again  where  water  would  flow)  of 
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Figure  21.  CN  Maximum  Numeric  Time  Oscillations 
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Figure  23.  CN  Maximum  Errors 
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least  errors,  though  remaining  parallel  for  both  cases,  became 
flatter  and  broader,  and  the  center  of  the  valley  shifted 
toward  higher  r  values.  Note  that  the  center  of  the  valley, 
as  well  as  the  10%  line,  are  roughly  parallel  and  linear.  By 
comparison  with  Figure  2  (top),  all  of  the  CN  FDMs  and  cases 
studied  satisfied  the  10%  criteria  of  the  full  matrix  method, 
which  predicted  the  linear  re  1  at i onship .  Finally,  the  range 
of  maximum  errors  decreases  with  the  increased  time  as 
expected,  since  the  transient  heat  equation  has  more  time  to 
reach  a  smaller  range  of  values  as  it  approaches  the 
temperature  of  the  boundary  conditions,  and  it  allows  the  FDM 
to  damp  the  solution  more  since  it  is  unconditionally  stable, 
and  approach  the  almost  steady  state  solution  with  no  discrete 
boundaries . 

Returning  to  the  standard  1-D  and  2-D  comperisons,  the  CN 
case  provides,  similar  to  the  EX,  clear  results  in  the  h-k 
domain.  For  time  oscillations,  very  high  r  correlations 
exist.  The  starting  time  oscillation  lines  start  cleanly  at 
about  r=2.0  and  r=1.0  for  the  1-D  and  2-D  cases  respectively, 
and  are  both  twice  as  high  as  predicted  by  the  coefficient 
method,  and  four  times  higher  than  predicted  by  the  more 
severe  full  matrix  predictions.  Hence,  the  time  oscillation 
predictions  are  both  conservative,  as  they  should  be,  but  by  a 
large  margin  as  well.  As  with  the  DF  spatial  oscillations, 
the  CN  spatial  oscillations  are  more  linear  where  they  start, 
and  do  not  correlate  well,  functionally,  with  the  predicted 
time  oscillations  based  on  stability  considerations. 
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Peaceman-Rachford  ADI  Results 


The  ADI  results  for  two  dimensions  are  show  in  Figures  24 
and  25.  Similar  to  the  DF  method,  the  ADI  method  exhibits 
complicated  behavior  in  the  h-k  domain.  In  order  to  correctly 
read  the  h-k  domain,  consideration  of  how  h  presented  on  the 
graph  relates  to  h  used  in  the  calculation  of  r  is  needed. 

The  ADI  program  breaks  the  user  given  time  step  into  two  time 
steps  to  perform  the  two  step  ADI  method,  where  each  of  the 
two  time  steps  has  an  r  that  is  one  half  that  which  the  user 
input  time  step  would  indicate.  Since  the  user  given  time 
step  is  graphed  on  the  vertical  axis,  the  observed  value  of  r 
on  the  h-k  domain  needs  to  be  divided  by  two  to  get  the  actual 
r  values  being  predicted. 

Note  that  unlike  the  other  methods  examined,  the  ADI  has 
a  valley,  with  high  r  correlation,  where  it  oscillates  at  a 
minimum  in  time,  with  shallow  hills  of  oscillation  on  either 
side.  The  shallow  hills  are  a  range  of  time  oscillations 
which  only  vary  from  1  to  4 ,  and  other  than  the  IM  method, 
show  the  smallest  range  for  any  of  the  1-D  or  2-D  cases 
examined  in  this  thesis.  So,  in  fact,  the  ADI  method  is  seen 
to  be  oscillatory  in  time  for  all  regions  in  the  h-k  domain 
for  this  sample  problem,  and  therefore  appears  to  be 
unconditionally  oscillatory  for  this  problem.  The  minimum 
valley  of  time  oscillation  was  found  to  occur  at  a  value  of 
one  time  oscillation,  located  between  r=0.75  to  r=2.0,  with 
less  than  two  oscillations  occurring  for  all  values  of  r<=2.0 
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Figure  25.  ADI  Maximum  Errors 


GAAAH  LBGCND 

HOA X IONTAL  AXIS: 
logl0(x-y  Atop) 
Length 

VIATICAL  AXIS: 
logl0(t  Atop) 
CoeondA 

X-AXIS : 

Total  Contour* 

92 

Start  1  ao  Lin 
1 

Start  1 ao  value 
9.73  58-02 
♦  /-  I  ao  oMaata 
2 . 000C-02 
X — Ax i A  Happing 
log  10 

Llsaar  Min  t  Hax 
2 . 4701400 
8 . 5271*01 
Linear  Avg  a  Sd 
1 . 2238*01 
4 . 1488*00 


in  the  domain.  The  coefficient  method  predicts  time 


e 


L 


■4 


oscillations  to  occur  at  r>l/2.  As  can  be  seen,  the 
coefficient  method  does  not  do  a  good  job  at  predicting  the 
complex  behavior  of  the  ADI  oscillations. 

Looking  at  the  spatial  oscillations,  the  more  familiar 
result  is  observed:  an  approximate  r  correlation,  though  more 
linear,  is  found.  Spatial  oscillations  start  at  a  line  that 
is  the  one  of  the  highest  in  the  h-k  domain  of  all  the  FDMs 
surveyed  in  this  study  (not  counting  no  lines  for  the  IM 
method).  The  starting  line  varies  from  above  r=3.0  to  r=6.0 
in  the  domain. 

Maximum  errors  of  the  ADI  method  never  become  larger  than 
85  degrees,  nor  smaller  than  2.5  degrees,  throughout  the 
domains  studied.  This  is  not  highly  accurate,  and  never  gets 
better  than  10%.  But  the  maximum  error  is  never  greater  than 
12.5  degrees  at  contour  line  51,  which  covers  most  of  the 
domain ,  and  includes  values  of  approximately  r<=4 . 0 .  Note 
that  the  DF  and  CN  methods  do  worse,  with  maximum  errors 
reaching  1200  and  260  degrees  in  the  domain,  respectively.  Of 
course,  the  DF  and  CN  methods  also  do  much  better,  with  the 
smallest  maximum  errors  reaching  2. 5x10-*  and  3. 2x10-^ 
degrees,  respectively.  And  since  12.5  degrees  is  the  range  of 
values  of  the  analytic  solution,  12.5  degrees  probably 
represents  unacceptable  error  in  the  final  solution,  as  it  can 
equal  100%  of  the  range  of  the  exact  solution.  Full  matrix 
predictions  were  not  completed  for  the  ADI  method,  so  they 
have  not  been  tested  against  the  observed  results. 
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VI.  Conclusions  and  Recommendations 


Cone lusions 

Based  on  the  results  of  this  study,  conclusions  will  now 
be  presented  concerning  two  topics:  1)  the  coefficient  and 

matrix  methods  of  predicting  stability  and  oscillations,  and 
2)  the  numerical  results  obtained  from  the  h-k  domain  graphs. 

The  Pade ’  method  was  found  to  be  equivalent  to  the  full 
matrix  method  for  predicting  stability,  where  only  the 
approach  to  deriving  the  amplification  matrix  differed. 
Further,  the  coefficient  method  p-ratio  was  found  to  be 
equivalent  to  the  expectation  value  of  the  matrix  method’s 
amplification  matrix  eigenvalues  for  two-level  FDMs .  In 
effect,  for  two-level  FDMs,  the  coefficient  method  p-ratios 
can  be  derived  from  the  matrix  method.  The  coefficient  method 
was  found  to  be  an  easy  to  use  and  apply  method  for  predicting 
stability  behavior,  though  not  as  powerful  as  the  matrix 
method.  However,  the  coefficient  method  should  do  well  at 
predicting  two-level  FDM  stability,  since  the  coefficient 
method  predictions  are  equivalent  to  the  expectation  value  of 
the  matrix  predictions. 

The  matrix  approach,  on  the  other  hand,  is  much  more 
difficult  to  apply,  especially  to  three-level  systems,  such  as 
the  DF.  Even  if  the  amplification  matrix  can  be  derived,  it 
may  be  very  complicated  to  analyze  to  determine  FDM  behavior. 


Further,  the  Pade ’  approach  to  deriving  the  amplification 
matrix  can  be  more  difficult,  if  not  impossible,  to  apply  to 
some  FDMs ,  such  as  the  DF ,  which  does  not  seem  to  be  easily 
factored  into  various  S,T  approx imant s .  Another  problem  with 
the  Pade’  approach  is  that  once  factored  into  several  S,T 
approximants  and  recombined,  such  as  in  the  case  of  the  ADI 
method,  one  has  to  interpret  how  to  replace  the  split  block 
matrices  with  eigenvalues  in  the  resulting  amplification 
matrix.  Even  when  the  eigenvalues  of  the  [B]  and  [C]  matrices 
are  known,  split  forms  involve  using  several  Pade’ 
approximants  that  do  not  reproduce  the  desired  FDM's  matrix 
unless  appropriate  truncations  are  carried  out.  In  general, 
it  is  simpler  and  safer  to  use  the  FDM  itself  to  construct  the 
amplification  matrix,  and  avoid  the  Pade’  approach,  when 
analyzing  a  given  FDM  for  stability. 

Due  to  the  more  mathematically  rigorous  nature  of  the 
matrix  method,  its  predictions,  though  severe,  should 
generally  be  more  exact  when  applied  to  an  infinite  set  of 
problems.  A  great  advantage  of  the  matrix  approach  is  that 
one  is  also  able  to  understand  and  predict  when  the 
oscillations  lead  to  errors  that  are  not  disastrous. 
Consequently,  a  less  severe  functional  relationship  in  the  h-k 
domain  can  be  derived,  since  in  practice,  small  errors  induced 
by  oscillations  can  still  lead  to  acceptable  solutions. 

The  h-k  domain  graphs  were  an  excellent  means  of 
comparing  predictions  based  on  stability  analysis  and  for 
finding  interesting  FDM  behavior.  The  h-k  domain  graphs 


showed  that  neither  the  coefficient  nor  matrix  method 
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predictions  were  significantly  violated  for  the  two-level  FDMs 
studied.  However,  unusual  DF  and  ADI  behavior  was  observed 
which  wasn't  predicted  by  the  coefficient  method.  Conclusions 
based  on  the  h-k  domain  graphs  will  be  covered  in  the 
following  order:  general  expected  and  unexpected  results, 
incorrect  diagnostics,  oscillatory  behavior  compared  to  the 
coefficient  and  matrix  method  predictions,  and  comments  on 
maximum  errors. 

Except  for  the  IM  method,  all  of  the  FDMs  surveyed  in  the 
h-k  domain  showed  varying  degrees  of  r  correlations,  as  was 
expected  from  the  major  role  r  plays  in  predicting  behavior. 
Also  expected  was  the  result  that  the  1-D  and  2-D  contour 
graphs ,  while  very  different  among  the  five  FDMs  examined, 
were  very  similar  for  the  same  FDM.  This  similarity  was 
maintained  regardless  of  the  different  problems,  final 
solution  times,  number  of  spatial  dimensions,  or  differences 
in  DF  jump  starts  used. 

Unexpected  results  were  observed  for  both  the  DF  and  ADI 
FDMs.  For  the  DF  FDM,  oscillations  were  not  found  in  regions 
with  large  values  of  r,  though  those  regions  had  large  errors. 
Moreover,  oscillations  were  found  in  regions  with  smaller  than 
predicted  values  of  r,  though  those  regions  had  errors  which 
were  quite  small.  Also,  the  DF  FDM  showed  unexpected  behavior 
when  the  analytic  solution,  IM,  and  CN  methods  were  used  to 
jump  start  the  first  time  level.  The  analytic  solution  jump 
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resulted  in  the  worst  time  oscillations  and  errors  in  the 


final  solution.  The  CN  jump  start,  an  FDM  known  to  be  able  to 
oscillate,  actually  resulted  in  i  he  least  time  oscillations, 
while  the  IM  jump  start  gave  the  most  accuracy.  Finally,  the 
ADI  method  always  displayed  time  oscillations  over  the  whole 
domain,  though  its  spatial  oscillations  were  more  normal  in 
appearance,  being  approximately  linear  and  only  occurring  high 
in  the  h-k  domain. 

Two  stability  diagnostics,  spatial  oscillations  and 
maximum  error  versus  time,  were  found  to  be  misleading  and 
incoriect.  Based  on  stability  analysis  of  FDMs ,  oscillations 
with  respect  to  time  (i.e.  time  steps)  are  predicted. 
Unfortunately,  spatial  oscillations  are  often  used  to  confirm 
stability  oscillation  predictions.  But,  in  this  study, 
spatial  oscillations,  except  for  the  EX  FDM,  were  found  to  be 
more  linear  (k=bh,  where  b  is  a  constant)  than  the  time 
oscillations  (k=rh2).  So  while  it  was  generally  observed  that 
spatial  oscillations  usually  occurred  in  regions  with  time 
osc i 1 1  at i ons ,  as  would  be  expected,  it  was  found  that  time 
oscillations  also  occurred  in  regions  of  the  domain  that  were 
not  oscillating  spatially.  Therefore,  spatial  oscillations 
used  to  prove  or  infer  stability  predictions  based  on  time 
oscillations  can  be  misleading.  At  best,  such  demonstrations 
often  require  that  larger  than  necessary  r  values  be  used  to 
find  oscillations,  and  thereby  tend  to  inflate  the  r  values 
that  stability  analysis  is  suppose  to  predict.  At  worst,  the 
use  of  spatial  oscillations  confuses  the  concepts  embodied  in 
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stability  predictions  as  well  as  providing  incorrect  results. 
Another  misleading  diagnostic  is  the  use  of  the  maximum  error 
versus  time.  Maximum  error  versus  time  is  not  a  good 
stability  (stable  decay  and  stable  oscillation)  indicator 
since  stability  is  a  necessary,  but  not  a  sufficient, 
condition  for  accuracy. 

Oscillatory  behavior  in  time  was  verified  to  exist,  as 
predicted  by  the  coefficient  method,  for  the  EX,  DF,  CN ,  and 
ADI  FDMs .  The  coefficient  method  did  not  predict  oscillations 
for  the  IM  method,  and  none  were  found.  Predicted  r  values 
for  the  start  of  oscillations  in  the  EX,  CN ,  and  DF  were  very 
good  at  predicting  values  below  which  oscillations  were  not 
occurring,  but  were  less  accurate  about  predicting  when  the 
oscillations  actually  began.  Such  conservative  predictions 
are  expected  as  other  test  problems  could  start  oscillating  in 
time  at  values  closer  to  the  predicted  r  values.  Overall,  the 
coefficient  method,  as  expected,  proved  to  be  good  at 
predicting  two-level  FDM  stability  behavior.  The  DF  and  ADI 
methods  exhibited  more  complex  behavior.  In  comparison,  the 
coefficient  method  only  showed  general  r  correlations  which 
only  approximately  predicted  the  oscillatory  behavior  observed 
in  the  DF  and  ADI  FDMs. 

The  matrix  method  was  completely  derived  and  applied  to 
the  IM,  EX  and  CN  methods.  As  with  the  coefficient  method, 
the  matrix  method  correctly  predicted  no  oscillations  for  the 
IM  method.  Since  the  predicted  r  values  (based  on  the  extreme 
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eigenvalues  of  the  matrix  method)  were  smaller,  usually  one 
half  the  value  of  the  predicted  coefficient  method  values,  the 
matrix  method  was  even  more  severe,  more  conservative  in  its 
predictions  than  the  coefficient  method.  This  severity  was 
expected,  since  from  the  derivations,  the  matrix  method 
stability  predictions  are  worst  case  predictions  that  should 
never  be  violated.  The  non-disastrous  oscillations 
predictions  proved  to  be  useful  in  determining  less  severe 
restrictions  on  h  and  k,  roughly  resulting  in  maximum  absolute 
errors  of  less  than  10%  of  the  final  range  of  temperatures  for 
the  problems  studied. 

Neither  the  coefficient  nor  the  matrix  method  was  proven 
to  be  substantially  incorrect  by  these  test  problems,  number 
of  dimensions,  or  final  solution  times  for  two-level  FDMs .  In 
addition,  the  coefficient  method  provided  acceptable  results 
without  the  severe  constraints  of  the  matrix  method.  Analysis 
of  the  DF  and  ADI  methods  using  the  matrix  method  could  not  be 

finished,  but  for  these  methods  the  eigenvalues  of  the 

amplification  matrices  indicate  that  more  complicated  behavior 
is  probable.  Therefore,  the  matrix  method  may  be  able  to 

explain  the  unexpected  results  seen  in  the  DF  and  ADI,  whereas 

the  simple  coefficient  method  could  not.  Further,  the 
eigenvalues  of  the  amplification  matrices  predicted  different 
behavior,  which  was  observed  in  the  test  problems,  for  the  1-D 
DF ,  2-D  CN ,  and  2-D  ADI,  whereas  the  coefficient  method  only 
predicted  the  same  behavior. 

The  worst  maximum  error  increased,  for  the  same  domains 
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and  for  both  the  1-D  and  2-D  cases,  in  the  following  order: 

IM,  ADI,  CN ,  DF ,  EX,  with  EX  having  the  worst  maximum  error 
due  to  its  conditional  stability  being  violated.  The  least 
maximum  error  generally  decreased  in  the  following  order: 

ADI,  IM,  CN,  EX,  DF ,  with  the  DF  have  the  least  maximum  error 
of  all  over  the  domain.  Generally,  the  CN  method  proved  to  be 
one  of  the  best  methods  for  small  maximum  errors  based  upon  a 
nearly  linear  restriction,  as  predicted  by  the  matrix  method, 
on  the  values  of  k/h,  versus  the  more  general  r  constraints 
based  on  k/h*.  For  h  and  k  of  any  size  in  the  domain,  the  IM 
method  provided  the  smallest  worst  maximum  error,  as  well  as 
having  the  smallest  range  of  maximum  errors  over  the  entire 
h-k  domain. 


Recommendat ions 

As  a  result  of  this  study,  the  following  are  proposed: 

1.  Finish  analyzing  the  1-D  DF ,  2-D  DF  and  ADI  FDMs 
using  the  full  matrix  method.  Explain  the  unexpected 
behavior,  not  predicted  by  the  coefficient  method,  observed  in 
the  DF  and  ADI  methods.  Relate  the  coefficient  method 
p-ratios  to  the  matrix  method  for  three-level  FDMs. 

2.  Confirm  or  reject  the  ADI  results  for  the  p-ratio, 
amplification  matrix  eigenvalues,  and  computer  graphs  in  this 
study.  Since  these  results  were  the  last  achieved  in  this 
study,  they  were  not  as  completely  tested  and  checked  as  oth-'r 
results,  and  represent  the  only  probable  areas  where 
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significant  errors  could  have  occurred  in  this  study. 

3.  Use  Gauss-Seidel  iteration  to  completely  jump  start 
the  DF  method  using  the  DF  method  itself,  which  is  what  the 
stability  predictions  are  based  on,  and  see  if  better  or  more 
predictable  results  are  obtained  then  by  using  the  analytic 
solution,  IM,  or  CN. 

4.  Compare  the  stability  predictions  based  on  the 
probabilistic  and  fourier  techniques  with  the  domain  graphs  of 
this  study. 

5.  Run  the  programs  for  other  problems  at  various  times 
to  see  if  the  less  rigorous  coefficient  method  can  be  found  to 
be  significantly  violated  in  its  predicted  non-osci 1 latory  r 
values.  Is  the  matrix  method  ever  violated? 

6.  Run  the  programs  to  see  if  the  probabilistic  method, 
as  applied  to  explicit  schemes,  is  ever  violated.  Also 
compare  the  probabilistic  method  with  the  matrix  method,  which 
considers  all  the  nodes,  to  determine  if  they  are  equivalent 
or  not.  Does  a  relationship  similar  to  that  of  the 
coefficient  and  matrix  methods  exist? 

7.  Develop  an  algorithm  to  plot  stability  in  the  h-k 
domain,  using  the  numerical  solution  being  bounded  as  a  trial 
diagnostic . 

8.  Apply  the  FDMs  studied  in  this  paper  to  problems  with 
derivative  or  mixed  boundary  conditions,  since  these  effect 
stability.  How  do  the  various  methods  of  predicting  stability 
compare  with  the  new  domain  results? 


9. 


Use  the  domain  graph  approach  to  explore  newly 


c 


developed  FDMs  to  see  if  they  oscillate  and  to  determine  if 
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h-k  Domain  Graphs  for  the  Slope  of  the 
Maximum  Error  Versus  Time 
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2-D  OuPort-Pranfcel ,  !■:  xl-yl-0,  *2-y2-l  and  tl,t2-0,.20s 

X-AXIS:  Slope  of  Haxiaiun  Error  Versus  Tiaiesteps 
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Main  Program  V  30 


Numerical  FDM  Solution 

t _ _ _ 

Exact  pde 
Solution 

i  i 

Decay/Growth 

i 

No  Orcillanon/Osdllation 

_ L_ 

Acceptable 
Error  Stats 

Important  1 

i 

i 

i - 

- 1 — 

Program  =»  uunQ 

Outputs 

uua(,) 

slpsd 

slpmax 

mtnosc 

mtnaosc 

snaosc 

(msnaosc) 

avg 

sd 

max  err 

Solid  Lines  —  1-D  and  2-D  paths  if  only  solid  lines  from  root 
Dotted  Lines  —  2-D  paths,  followed  instead  of  solid  lines  from  a  root 
[  ]  —  optional,  if  required 
Underline  —  function 
No  Underline  —  procedure 
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Definitions  of  Structure  Diagram  Terms 


fdmuun: 

pdeuua: 

trenderr: 

maxtimosc: 

spaosc: 

errorstats: 

tincc: 

bcxl: 

bcx2: 

[tridi] 

exactfunc: 

errorstats: 

Unfit: 

timosc: 


calculates  and  stores  all  Finite  Difference  Method  (FDM)  solutions  in 
uun  (ix,  it)  for  all  nodes,  all  timesteps 

calculates  and  stores  all  partial  differential  equation  (pde)  solutions  in 
uua  (ix,  it)  for  all  nodes,  all  timesteps 

calculates  sd  and  max  error  slopes  of  absolute  (numeric  solution  - 
analytic  solution)  errors  at  each  timestep  versus  that  timestep 
calculates  maximum  N  and  N  -  A*  time  oscillations  from  all  spatial 
nodes 

calculates  total  spatial  oscillations  at  any  timestep  —  here  it  is  used  for 
final  timestep,  nt 

calculates  error  statistics  tor  all  nodes  (includes  BCs)  for  anv 
timestep  —  here  it  is  used  for  final  timestep,  nt 

temperature  initial  condition  (IC) 
left  boundary  condition  (BC) 
right  boundary  condition  (BC) 
tridiagonal  solver  (if  needed) 

exact  analytic  series  solution  to  pde 

calculates  error  statistics  (aug,  sd,  max  absolute  errors)  over  all  nodes 
for  anv  timestep 

least  squares  curve  fitting  routine  to  calculate  error  slopes  (sd  &  max) 
for  trending  error  behavior 

calculates  total  time  oscillations  at  any  spatial  node  for  N  and 
N  -  A*  solutions 


*  N  —  Numeric  Solution,  A  —  Exact  Analytic  Solution 


Calculation  of  Program  Run  Time 


-(x2-x0 


Qa-'Q 


Approximate  total  node  points  solved  in  h-k  domain: 


1D.  (max  n^)  (maxnj 
4 


3  2 

(max  (maxn!) 

2-D: - - -  for  Ax  =  Ay  plane  (used  in  this  thesis,  with  nx  =  ny  also) 


^  n.  (maxn0  (maxny)  (max  n  j 
8 


for  all  Ax  and  Ay 


Note:  ICC  can  solve  for  1.125  million  node  points  per  hour  on  average,  with 
0.675  to  2.25  million  node  points  per  hour  for  ±  3  standard  deviations. 


Example:  1-D:  nx  =  3  tolOO,  nt  =  3  to  100  =>*  25  million  points  *  22  hours  to  run 
2-D:  nx  =  ny  =  3  to  25,  nt  =  3  to  100  =>  =  26  million  nodes  ~  23  hours 


Appendix  C: 


3ilitv  DenvatK 


1-D  EX  Example 


End  =  0‘ 


oc  =  1  =  Thermal  diffusmty 
Rod  =  1* 


End  =  O' 


T  =  M  At  U(N,n) 

U(0,n  +  1)  I 


L  =  N  Ax 


~W~a  .  2 

dx 

continuous 

dU  U  (x,  t  +  At)  -  U  (x,  t) 

^  At 

to 

l)(x  +  Ax,t)-2U(x,t)  +  u(x-Ax,t) 

2  2 

discrete 

d*  Ax 

Substitute  derivative  approximations  into  pde  and  rearrange  terms  to  get  explicit  form: 

U  (x,  t  +  At)  =  rU  (x  -  Ax,  t)  +  (1  -  2r)  U  (x,  t)  +  rU  (x  -  Ax,  t) 

which  can  be  written  in  matrix  form  as: 


™l-2r  r 
r  l-2r  r 

"  Ul.n  ' 

u2.„ 

- 1 

=  £ 
ss 

_ 1 

'  U1<n+1  ’ 
U2.n+1 

- 

- 

+ 

- 

= 

— 

r  l-2r  r 
r  1— 2r 

U  N-2,n 
UN-l.n 

_ 1 

U  N-2.H+1 
Un-i.„+i 

Coefficient  Method: 


P  = 


U  (x,  t  +  At) 
U  (x,  t) 


1  -2r 
1 


=  1  -  2r 


r  =  0  to  ^/2  stable  decay 

T~Vl to  1  stable  osc 
r>  1  unstable  osc 


Matrix  Method: 


Tridi  Matrix  Zs  =  b  +  2  Vac  cos s  =  1,  2, ...  N-l 


S7t 


.  2  SJt 


|is = Gs  =  1  -  2r  +  2r  cos  =  1  —  4r  sin  ^ 

=  1  -  2r  =  coefficient  method  p-ratio 
Extremes:  M-j  ~  1,  HN_  j  » 1  -  4r  for  large  N 


r  =  0  to  I/4  stable  decay  | 

1  —  4r  /  r=y^to]/2  stableosc  \ 
^/2  unstable  osc  I 


r> 


Non-disastrous  oscillations  condition: 


for  eigenvalue  extremes 


M-n-  1 


r  < 


or  r  <  approximately  for  large  N 


Matrix  Non-Disastrous  Oscillations 


Condition: 


-Hj  <Mn- 1  <^, 


.  ,  .in  .  2N-I  7C  2  ic 

-  1  +  4r  sin  <  1  -  4r  sin  1 — ==-= —  <  1  -  4r  sin 

2N  2N  2N 


RHS:  sin 


.  2  (N  —  1}  7C  .  2  n 


2N 


>  sin  2^-  is  always  satisfied 


,  .  2  7t  „  .  .  2(N  -  lht 

LHS:  4r sin  2^-<2-4rsin  — ^ — 


LET: 


.  2  n  7t 
sin  2N~T^2 


4N  \  for  large  N 


.  2(N-1)k 
sin 


2N 


7t 


THEN:  4r - <2-4r 

4N2 


iL  +  4 

\N2  / 


<2 


r< 


iL  +  4 

N2 


therefore 


r  <■ 


2  + 


7C 


2Ni 


condition  for  non -disastrous  oscillations 


NOTE:  Approaches  r  <  1/2  for  large  N,  where  (ls  varies  from  -1  to  1  for  (iN_  j  =  1 
Therefore,  agrees  with  the  condition  based  on  I  eigenvalues  I  <  1. 


Derivation  of  the  Two  Dimensional  EX  Full  Matrix 
Parameters  using  the  Pade'  Method 


Let  mnj,  *  Gm  n  *  Amplification  Matrix  Eigenvalues 

From  the  2-D  EX  Pade'  Appro ximant:  4^=  1  +  kZm>n 
„  4  [  .  imit  .  2  n«  \ 

where  Z^-  —  ^  ^j  +  sin 


and m  =  1  (1) N -  1;  n  =  l(l)N-l 


,  .  [  .  imn  .  2ntt  ] 

2N +sm  2N ) 


,  then 


=  1  -  4r  =  coefficient  method  p-ratio 
Extremes:  lAu*  1,  1  -  8r 


1  —  8r 


stable  decay\ 
stable  osc 
unstable  osc 


) 


for  eigenvalue  extremes 


For  non-disastrous  oscillations,  the  high  frequency  components 
must  decay  faster  than  the  low  frequency  components  so  that 


Condition:  Hn-i.N-i  <  Hi.i 


RHS:  is  always  satisfied 


LHS:  —  1  -  kZjj<  1  +  j 

"kZi  i<2  +kZN_i  n_i 


LET: 


Zll=-A[  JL_+JL]  \ 

h  \^4N  4N  )  V  for  large  N 

^N-l,N-l~  -  ~2  J 


THEN:  4r — -<2-8r 
4N2 


r  ^-  +  8  <2 
N2 


r< — - therefore 

271 

”  +  8 
N 


condition  for  non-disastrous  oscillations 


4  +  — 
N2 


e 


l-DIM 


1-DCN 


1-DDF 


2-DDF 


2-DADI 


Appendix  D: 


Amplification  Matrix  Eigenvalues 


9  S7t 

l-DEX  l-4rsin 

2N 


2-DEX 


,  A  .  2  171  2  J7t 

l-4r  sm  ™ 


2-DIM 


i  ^  A  ■  2  SK 

1  4rSU1  2N 


,  ,  .  2  17t  .  2  J7C 

1  +  4r  sin  ^i-r  +  sin  ~ 
2N  2N 


i  o  •  2  S7t 

1  -  2r  sm 


2-DCN 


,  „  .  2  17C  .  2  JJt 

1_2r  Sln  2N+SU1  2N 


,  „  .  2  sn 

1  +  2r  sin  ^ 


,  ~  .  2  ITt  .  2  J7t 

l+2r  sm2N+s“  2N 


o  S7C  /  ,  2  .  2  S7C 

2r cos ± 'Y  l-4r  sin  -jj- 


1  +  2r 


_  17t  lit  + 

2r  cos  -rr  +  cos  Vr  - 
N  N/ 


,  A  2  a  171  J71 

l_4r  4  —  cos  +  cos  vr 
N  N 


1  +  4r 


,  -  .  2J7t  ,  .  .  11K 

1  -  2r  sin  2N  1  ~  2r  Sin  2N 


,  -  •  2  lit  ,  _  .  2  J7t 

1  +  2r  sm  ^  1  +  2r  sin 


where  s  =  1,  2, ...  N  -  1 
i  =  1,  2, ...  N  -  1 
j  =  1, 2, ...  N  -  1 


*  i  * 


Expectation  Values 


/.  2  S7l\  /.  2S7c\  /  2  Slt\  /  2Slt\_l 

\SU1  N/  =  \sin  2N/  =  \COS  7T/  =  \COS  2N/  '  2 


cos^)  =  0 


— /  =  \cos2|^-)  =  2  N  >  2  s  =  1, 2, ...  N  - 1 


N  >  2  i  =  1, 2, ...  N  -  1 


.  2  nt 
sin  XTT  +  sin 
v  2N 


.  2  jit  \  /  .  2  ill  \  ,  /  .  2  jit  \  , 

m  2N/  =  V”  2N/  Vm  2N/ 


lit  Jit 

|cos^  +  cosJ— I 


=  1  for  large  N 
=  0  N  =  2 


N£2  i  =  1,  2, ...  N  -  1 
j  =  1, 2, ...  N  -  1 


i=  1,2,  ...N  -  1 
j  =  1,  2, ...  N  -  1 


(  )  indicates  expectation  value  of  enclosed  quantity,  defined  as: 


,  N-i 

<kzC>*Njrr 

8=  1 


or 


f- 


<kZ,)- 


4r  sin2  0  d8 


/. 


% 


for  large  N 


d0 
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Comparison  of  Coefficient  P-Ratio  and 
Expectation  Values  of  Amplification  Matrix  Eigenvalues 


FDM 

Coef  P-Ratio 

(Gs)  or  (Gm>n) 

1-DEX 

1  -2r 

Same 

2-DEX 

1  -4r 

Same 

1 

1-DIM 

1  +  2r 

Same 

1 

2-DIM 

1  +  4r 

Same 

1  -r 

1-DCN 

1  +  r 

Same 

1  -2r 

*  2-D  CN 

1  +2r 

Same 

1  -2r 

±  V  1  -  2r2 

*  1-D  DF 

1  +  2r 

1  +  2r 

1  -4r 

±  V  1  -  12r 

**  2-D  DF 

1  +  4r 

1  +4r 

1  -  2r 

1  -2r  +  r2 

*  2-D  ADI 

1  +2r 

1  +  2r  +  r 

* 

Have  same  Coef  P-Ratio 

** 

±  V  1  -  16r2 

f or  N  =  2 

1  +4r 
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